{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multi-point stress approximation (MPSA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Porepy supports mpsa discretization for linear elasticity problem:\n",
    "\\begin{equation}\n",
    "\\nabla\\cdot \\sigma = -\\vec f,\\quad \\vec x \\in \\Omega\n",
    "\\end{equation}\n",
    "where $\\vec f$ is a body force, and the stress $\\sigma$ is given as a linear function of the displacement\n",
    "\\begin{equation}\n",
    "\\sigma = C:\\vec u.\n",
    "\\end{equation}\n",
    "\n",
    "The convention in porepy is that tension is positive. This means that the Cartesian component of the traction $\\vec T = \\sigma \\cdot \\vec n$, for a direction $\\vec r$ is positive number if the inner product $\\vec T\\cdot \\vec r$ is positive. The displacements will give the difference between the initial state of the rock and the deformed state. If we consider a point in its initial state $\\vec x \\in \\Omega$ and let $\\vec x^* \\in \\Omega$ be the same point in the deformed state, to be consistent with the convention we used for traction, the displacements are given by $\\vec u = \\vec x^* - \\vec x$, that is, $u$ points from the initial state to the finial state.\n",
    "\n",
    "To close the system we also need to define a set of boundary conditions. Here we have three posibilities, Neumman conditions, Dirichlet conditions or Robin conditions, and we divide the boundary into three disjont sets $\\Gamma_N$, $\\Gamma_D$ and $\\Gamma_R$ for the three different types of boundary conditions\n",
    "\\begin{equation}\n",
    "\\vec u = g_D, \\quad \\vec x \\in \\Gamma_D\\\\\n",
    "\\sigma\\cdot n = g_N, \\quad \\vec x \\in \\Gamma_N\\\\\n",
    "\\sigma\\cdot n + W \\vec u= g_R,\\quad \\vec x\\in \\Gamma_R\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "To solve this system we first have to create the grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import porepy as pp\n",
    "\n",
    "# Create grid\n",
    "n = 5\n",
    "g = pp.CartGrid([n,n])\n",
    "g.compute_geometry()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to define the stress tensor $C$. In porepy the constutitive law,\n",
    "\\begin{equation}\n",
    "\\sigma = C:u = 2  \\mu  \\epsilon +\\lambda  \\text{trace}(\\epsilon) I, \\quad \\epsilon = \\frac{1}{2}(\\nabla u + (\\nabla u)^\\top)\n",
    "\\end{equation}\n",
    "is implemented, and to get the tensor for this law we call:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create stiffness matrix\n",
    "lam = np.ones(g.num_cells)\n",
    "mu = np.ones(g.num_cells)\n",
    "C = pp.FourthOrderTensor(mu, lam)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we need to define boundary conditions. We set the bottom boundary as a Dirichlet boundary, and the other boundaries are set to Neuman."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define boundary type\n",
    "dirich = np.ravel(np.argwhere(g.face_centers[1] < 1e-10))\n",
    "bound = pp.BoundaryConditionVectorial(g, dirich, ['dir']*dirich.size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We discretize the stresses by using the multi-point stress approximation (for details, please see: E. Keilegavlen and J. M. Nordbotten. “Finite volume methods for elasticity with weak symmetry”. In: International Journal for Numerical Methods in Engineering (2017))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now define the values we put on the boundaries. We clamp the bottom boundary, and push down by a constant force on the top boundary. Note that the value of the Neumann condition given on a face $\\pi$ is the integrated traction $\\int_\\pi g_N d\\vec x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_faces = np.ravel(np.argwhere(g.face_centers[1] > n - 1e-10))\n",
    "bot_faces = np.ravel(np.argwhere(g.face_centers[1] < 1e-10))\n",
    "\n",
    "u_b = np.zeros((g.dim, g.num_faces))\n",
    "u_b[1, top_faces] = -1 * g.face_areas[top_faces]\n",
    "u_b[:, bot_faces] = 0\n",
    "\n",
    "u_b = u_b.ravel('F')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We discretize this system using the Mpsa class. We assume zero body forces $f=0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of ncasym eliminated:  8\n",
      "max igrad:  2.0\n"
     ]
    }
   ],
   "source": [
    "mpsa_class = pp.Mpsa(\"mechanics\")\n",
    "f = np.zeros(g.dim * g.num_cells)\n",
    "\n",
    "specified_parameters = {\"fourth_order_tensor\": C, \"source\": f, \"bc\": bound, \"bc_values\": u_b}\n",
    "data = pp.initialize_default_data(g, {}, \"mechanics\", specified_parameters)\n",
    "A, b = mpsa_class.assemble_matrix_rhs(g, data)\n",
    "\n",
    "u_class = np.linalg.solve(A.A, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can plott the y_displacement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/rbe051/anaconda3/envs/porepy/lib/python3.6/site-packages/mpl_toolkits/mplot3d/axes3d.py:738: UserWarning: Attempting to set identical bottom==top results\n",
      "in singular transformations; automatically expanding.\n",
      "bottom=0.0, top=0.0\n",
      "  'bottom=%s, top=%s') % (bottom, top))\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxkAAAKhCAYAAAA8DWNTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xl8leWd///3lXNCICEBWQ0hsgiy\nrwmLS60LKLUtFKWIisIgjp0yM3X0O6MtSml11Nqx1Rlt52er0ypfl7bfig4qlcXWqiCEbUQQEQgQ\nEgLIluRwkrNcvz+yDCHnJHBzk/sOeT0fj/PgnHPfufPJ7THcH97XdV/GWisAAAAAcEuK1wUAAAAA\nOL/QZAAAAABwFU0GAAAAAFfRZAAAAABwFU0GAAAAAFfRZAAAAABwFU0GAAAAAFfRZAAAAABwFU0G\nAAAAAFcFz3B/lgcHAADAuWa8LuB09DPGhjyuoUT6k7V2ksdlNHCmTQYAAAAASSFJd3tcw0Kpi8cl\nJMRwKQAAAACuIskAAAAAHDDiYjoZkgwAAAAArqLJAAAAAOAqEh4AAADAASMp1esifIokAwAAAICr\nSDIAAAAAB5j4nRxJBgAAAABX0WQAAAAAcBUJDwAAAOAAE7+TI8kAAAAAWiFjTCdjzDJjzPaaPy9I\nsM9IY8wqY8ynxpj/McbcfDrHpskAAAAAHKid+O3l4yw9IGmFtba/pBU1r08VknSHtXaIpEmSnjLG\ndGzqwDQZAAAAQOs0RdJva57/VtK3Tt3BWvu5tXZ7zfNiSQckdW3qwMzJAAAAAFquLsaYgpNeP2et\nfe40v7a7tbak5vl+Sd0b29kYM1ZSG0k7mjowTQYAAADggE8mfh+y1uYn22iMWS7pwgSb5p/8wlpr\njTG2keNkS3pJ0ixrbbypomgyAAAAgPOUtXZCsm3GmFJjTLa1tqSmiTiQZL8sSW9Jmm+tXX0635cm\nAwAAAHDgPFjx+01JsyQ9XvPnG6fuYIxpI+l1SS9aa/9wugdm4jcAAADQOj0uaaIxZrukCTWvZYzJ\nN8b8umaf6ZKulDTbGLOx5jGyqQMba5MOvUrkjHYGAAAAHDBeF3A6+hhjf+hxDX8jrWtsToZXWnjC\nAwAAAHjDJxO/fYnhUgAAAABcRZIBAAAAOHAeTPw+Z0gyAAAAALiKJgMAAACAq0h4AAAAAAeY+J0c\nSQYAAAAAV9FkAAAAAHAVw6UAAAAABxgulRxJBgAAAABXkWQAAAAADnExnRhJBgAAAABX0WQAAAAA\ncBUJDwAAAOAAE7+TI8kAAAAA4CqSDAAAAMABIy6mkyHJANCitG/fXjt37ky47Te/+Y2uuOKKZq4I\nAACciiYDQAMvv/yy8vPz1b59e2VnZ+trX/uaPvjgA8fHM8boiy++qPdeWVmZ7r33XvXu3VsZGRm6\n6KKLNG3aNH388ceNHqu8vFx9+/Z1VEdVVZUWLlyo/v37KyMjQ71799acOXNUWFjo6HjnysKFCzVz\n5kyvyzinCgsLZYxRNBr1uhQAwDlAkwGgnp/97Ge655579IMf/EClpaXas2ePvvvd7+qNN94442Ml\nu4CsrKzUNddco08++URLlizR8ePHtXXrVs2YMUPvvPPOGR3rTEybNk1vvvmmXn75ZR07dkybNm1S\nXl6eVqxYcdbHhvtoQAD4Xe3Eby8ffkWTAaDOsWPHtGDBAj377LO68cYblZGRodTUVH3zm9/UT3/6\nU0nSmjVrdOmll6pjx47Kzs7W3//936uqqqruGMYYPfvss+rfv7/69++vK6+8UpI0YsQItW/fXq+9\n9ppeeuklFRUVafHixRo6dKgCgYAyMjI0bdo0LVy4MOmxat+rTUW+/PJLTZ48WVlZWRo7dqx27NiR\n9Gdbvny5li1bpjfeeENjxoxRMBhUhw4dNG/ePN15552SpOLiYk2ePFmdOnVSv3799Ktf/aru6xcu\nXKhvf/vbmjlzpjIzMzVs2DB9/vnneuyxx9StWzfl5ubq3Xffrdv/qquu0ve//32NHTtWWVlZmjJl\nig4fPixJ+vOf/6yePXvWq693795avny5li5dqkcffVSvvfaa2rdvrxEjRtT9t7nzzjuVnZ2tnJwc\nPfjgg4rFYgl/1hMnTmjWrFm64IILNGjQID3xxBP1vl9xcbFuuukmde3aVX369NG///u/122rrKzU\nPffcox49eqhHjx665557VFlZWa/uJ554Qt26dVN2drYWL16st99+W5dccok6deqkRx99tO5Y8Xhc\njz/+uC6++GJ17txZ06dPrzsHtZ+Ljh07qn379lq1apV+85vf6PLLL9c//dM/qXPnzlqwYIE6deqk\nTz75pO6YBw4cUHp6ug4ePJj0vzUAwHs0GQDqrFq1SuFwWFOnTk26TyAQ0M9//nMdOnRIq1at0ooV\nK/SLX/yi3j6LFy/Wxx9/rC1btuj999+XJG3atEnl5eW6+eabtXz5cl1//fXKyMhosqaTj3WqefPm\nqW3btiopKdELL7ygF154Ielxli9frrFjxyo3NzfpPjNmzFDPnj1VXFysP/zhD/rBD36glStX1m3/\n7//+b91+++06cuSIRo0apeuvv17xeFz79u3TggULdPfdd9c73osvvqgXXnhBJSUlCgaD+sd//Mcm\nf95JkybpBz/4gW6++WaVl5dr06ZNkqTZs2crGAzqiy++0IYNG/Tuu+/q17/+dcJj/OhHP1JhYaF2\n7typZcuWadGiRXXb4vG4vvnNb2rEiBHat2+fVqxYoaeeekp/+tOfJEn/+q//qtWrV2vjxo3atGmT\n1qxZo0ceeaTu6/fv369wOKx9+/bpxz/+se666y4tWrRI69at01//+lc9/PDD2rVrlyTpP/7jP7R4\n8WL95S9/UXFxsS644ALNmzdPkuo+F0ePHlV5ebkuvfRSSdLHH3+svn37qrS0VA899JBmzJhRr/5X\nXnlF1157rbp27drkuQSAc6124reXD9+y1p7JA8B5bNGiRbZ79+5n9DU///nP7be+9a2615LsihUr\n6u0jyW7fvr3u9bXXXmvvv//+utcbNmywHTp0sJmZmfaSSy45rWNFo1EbDAbt1q1b67Z9//vft5df\nfnnCOufOnWtvvvnmpD/Hnj17bEpKij1+/Hjdew888ICdNWuWtdbaH/7wh3bChAl12958802bkZFh\no9Gotdba48ePW0n2yJEj1lprv/rVr9b7GT/99FObmppqo9Gofe+992xOTk6979+rVy+7bNmyuu91\n22231W3bv3+/bdOmjQ2FQnXvvfzyy/aqq65K+LP06dPHLl26tO71r371q7rvt3r1apubm1tv/0cf\nfdTOnj3bWmtt37597VtvvVW3benSpbZXr17WWmvfe+8927Zt2wY/8+rVq+v2Hz16tH399dettdYO\nHDjQLl++vG5bcXGxDQaDNhKJ2F27dllJNhKJ1G3/r//6rwa11dYbj8ettdbm5eXZ1157LeHPDeC8\ncqbXqJ48Bkp2tccPSQVen4dED183QACaV+fOnXXo0CFFo1EFg4l/PXz++ee69957VVBQoFAopGg0\nqry8vHr7NJYW1H6fkpKSutcjR47U0aNHtXz5cs2dO/e0jnXw4EFFo9F623v16tXo9/z888+Tbi8u\nLlanTp2UmZlZ73gFBQV1r7t37173vF27durSpYsCgUDda6l6YnrHjh0b1N6rVy9FIhEdOnQoaQ3J\n7N69W5FIRNnZ2XXvxePxpOemuLi43raTn+/evVvFxcV1NUpSLBbTV77ylbqvPfk89urVS8XFxXWv\nO3fu3OBnPvW8lJeX132vqVOnKiXlf0PzQCCg0tLSpD/rqT/TuHHjlJ6erj//+c/Kzs7WF198ocmT\nJyf9egCAPzBcCkCdSy+9VGlpaVq8eHHSff7u7/5OAwcO1Pbt23X8+HE9+uijstbW28cY0+j3ufba\na/Xuu++qoqKiyZqSHatr164KBoPau3dv3Xt79uxJepwJEyZozZo1KioqSri9R48eOnz4sMrKyuod\nLycnp8kakzm1ttTUVHXp0kUZGRkKhUJ122KxWL05Bqf+zLm5uUpLS9OhQ4d09OhRHT16VMePH9en\nn36a8PtmZ2fX+zlPriM3N1d9+vSpO87Ro0dVVlamt99+W1L1edi9e3e9unv06OHo58/NzdU777xT\n73uFw2Hl5OQk/e+a6P1Zs2Zp0aJFeumllzRt2jS1bdvWUT0A4DYmfidHkwGgTocOHfTjH/9Y8+bN\n0+LFixUKhRSJRPTOO+/oX/7lXyRV33o2KytL7du312effaZf/vKXTR63e/fu9da2uOOOO5Sdna2p\nU6dq8+bNisViCofD9VKDpgQCAd14441auHChQqGQtmzZot/+9rdJ958wYYImTpyoqVOnat26dYpG\noyorK9N//ud/6oUXXlBubq4uu+wyff/731c4HNb//M//6Pnnnz+rW8kuWrRIW7ZsUSgU0oIFCzRt\n2jQFAgFdcsklCofDeuuttxSJRPTII4/UTa6Wqs9XYWGh4vG4pOqm4brrrtN9992n48ePKx6Pa8eO\nHfrLX/6S8PtOnz5djz32mI4cOaJ9+/bpmWeeqds2duxYZWZm6ic/+YlOnDihWCymzZs3a+3atZKk\nW265RY888ogOHjyoQ4cO6cc//rHjc/Cd73xH8+fPr2taDh48WHeXsq5duyolJSXpmicnmzlzpl5/\n/XUtWrRId9xxh6NaAADNiyYDQD333Xeffvazn+mRRx5R165dlZubq2eeeUbf+ta3JEn/9m//ppdf\nflmZmZm66667dPPNNzd5zIULF2rWrFnq2LGjfve736lt27Z67733NHjwYH39619XVlaWBgwYoLVr\n1+p3v/vdadf6zDPPqLy8XBdeeKFmz56tv/mbv2l0/z/84Q+64YYbdPPNN6tDhw4aOnSoCgoKNGHC\nBEnVk4oLCwvVo0cPTZ06VT/60Y/qtjlx++23a/bs2brwwgsVDofr7uLUoUMH/eIXv9DcuXOVk5Oj\njIyMend/+va3vy2pemjS6NGjJVVPIq+qqtLgwYN1wQUXaNq0afWGnJ1swYIF6tmzp/r06aMJEyZo\n2rRpSktLk1TdnC1ZskQbN25Unz591KVLF82dO1fHjh2TJD344IPKz8/X8OHDNWzYMI0ePVoPPvig\no5//e9/7niZPnqzrrrtOmZmZGj9+fN06KOnp6Zo/f74uv/xydezYUatXr056nNzcXI0ePVrGmLph\nXQDgB0z8Ts6cOsyhCWe0MwC0VldddZVmzpzZYI6JF375y1/q1VdfTZp8tARz5sxRjx496t3pCsB5\nrfFxtz4xxBj7isc1jJDWWWvzPS6jAZIMADjPlJSU6MMPP1Q8Hte2bdv05JNPNnpbYr8rLCzUH//4\nx7r1TAAA/keTAQDnmaqqKt19993KzMzUNddcoylTpui73/2u12U58tBDD2no0KH653/+Z/Xp08fr\ncgCgHiZ+J8dwKQAAAPhNixguNdQY+3uPaxjs0+FSfp4vAgAAAPhW7cRvNMRwKQAAAACuoskAAAAA\n4CoSHgAAAMCB2onfaIgkAwAAAICraDIAAAAAuIrhUgAAAIADDJdKjiQDAAAAgKtIMgAAAACHuJhO\njCQDAAAAgKtoMgAAAAC4ioQHAAAAcMBISvX6ajrq8fdPgiQDAAAAgKu87r0AAACAFskYKej11TRJ\nBgAAAIDWgCYDAAAAgKu8DngAAACAFskYKTXgdRX+RJIBAAAAwFUkGQAAAIADvpj47VMkGQAAAABc\nRZMBAAAAwFUEPAAAAIADvljx26dIMgAAAAC4it4LAAAAcMJI4ha2CZFkAAAAAHAVTQYAAAAAVzFc\nCgAAAHDCiKvpJEgyAAAAALiK3gsAAABwgiQjKZIMAAAAAK6iyQAAAADgKgIeAAAAwCmuphMiyQAA\nAADgKpoMAAAAAK4i4AEAAACcMJICXhfhTyQZAAAAAFxFkgEAAAA4wToZSZFkAAAAAHAVTQYAAAAA\nVxHwAAAAAE4wXCopkgwAAAAArqL3AgAAAJziFrYJkWQAAAAAcBVNBgAAAABXMVwKAAAAcIKJ30mR\nZAAAAABwFb0XAAAA4ARJRlIkGQAAAEArZIzpZIxZZozZXvPnBY3sm2WMKTLGPHM6x6bJAAAAAFqn\nByStsNb2l7Si5nUyD0t6/3QPTJMBAAAAOBXw+HF2pkj6bc3z30r6VqKdjDF5krpLevd0D0yTAQAA\nALRO3a21JTXP96u6kajHGJMi6UlJ/+dMDsxUFQAAAMAJf0z87mKMKTjp9XPW2udqXxhjlku6MMHX\nzT/5hbXWGmNsgv2+K+lta22RMea0i/L+tAAAAABw6pC1Nj/ZRmvthGTbjDGlxphsa22JMSZb0oEE\nu10q6SvGmO9Kai+pjTGm3Frb2PwNmgwAAACglXpT0ixJj9f8+capO1hrb6t9boyZLSm/qQZDoskA\nAAAAnPHHcKmz8bik3xlj7pS0W9J0STLG5Ev6jrV2rtMDG2sTDb1K6ox2BgAAABw4/cH/HsrPMrZg\nrLc1mBVa19hwKa+07N4LAAAA8ErLTzLOGW5hCwAAAMBVNBkAAAAAXEXAAwAAADh19qtun5dIMgAA\nAAC4iiYDAAAAgKsYLgUAAAA4wd2lkiLJAAAAAOAqei8AAADACZKMpEgyAAAAALiKJgMAPLJ27VoN\nHz5c4XBYFRUVGjJkiDZv3ux1WQAAnDUCHgDwyJgxYzR58mQ9+OCDOnHihGbOnKmhQ4d6XRYA4HQZ\nsU5GEsZaeyb7n9HOAIDGVVVVacyYMWrbtq0++ugjBQL8bQUAqr589738TsYWTPC2BvN7rbPW5ntb\nRUMkGQDgoS+//FLl5eWKRCIKh8PKyMjwuiQAwOli4ndSzMkAAA/dfffdevjhh3Xbbbfp/vvv97oc\nAABcQe8FAB558cUXlZqaqltvvVWxWEyXXXaZVq5cqWuuucbr0gAAOCvMyUCziUajMsYw5hwAADSl\nZczJ6Gxswde9rcG85M85GQyXQrOJRqPauXOnotGo16UAAADgHGK4FJpVUVGRevbsKUkKBvn4AQCA\nFoxb2CZFkoFmZ4xRJBIh0QAAADhP0WSg2RljZIzRrl27aDQAAADOQ4xXgSeMMQydAgAALRvrZCRF\nkgFPkWgAAACcf+i94KlTE41AICBjWsRd6wAAQGtHkpEUSQaaTUlJiSorKxNuq50MvnPnTp3h2i0A\nAADwGZoMNJvu3bsrFotp+/btDRqJ2snge/bsUTQaPetGY+nSpRowYID69eunxx9//KyOdT6YM2eO\nunXrpqFDh3pdii/s3btXV199tQYPHqwhQ4bo6aef9rokT4XDYY0dO1YjRozQkCFD9MMf/tDrknwj\nFotp1KhR+sY3vuF1Kb7Qu3dvDRs2TCNHjlR+vu/W/gLgIzQZaDYpKSlKT09XVVWVtm7d2mB77TCp\naDR6Vo1GLBbTvHnz9M4772jLli165ZVXtGXLlrOqvaWbPXu2li5d6nUZvhEMBvXkk09qy5YtWr16\ntZ599tlW/RlJS0vTypUrtWnTJm3cuFFLly7V6tWrvS7LF55++mkNGjTI6zJ85b333tPGjRtVUFDg\ndSmAPwQ9fvgUTQaa3eDBgxUIBHTixAnF4/EG240xKiwsdNxorFmzRv369VPfvn3Vpk0bzZgxQ2+8\n8YYbpbdYV155pTp16uR1Gb6RnZ2t0aNHS5IyMzM1aNAg7du3z+OqvGOMUfv27SVJkUhEkUiEuVGq\nXjz0rbfe0ty5c70uBQBaHJoMNDtjjAYMGKCUlBRt3Lgx4dCpoqIix4nGvn37lJubW/e6Z8+erfoC\nEo0rLCzUhg0bNG7cOK9L8VQsFtPIkSPVrVs3TZw4sdWfD0m655579MQTTyglhb8qaxljdN111ykv\nL0/PPfec1+UA8DF+c8IzaWlp6tKli0KhUMJb2J5togE0pby8XDfddJOeeuopZWVleV2OpwKBgDZu\n3KiioiKtWbNGmzdv9rokTy1ZskTdunVTXl6e16X4ygcffKD169frnXfe0bPPPqv333/f65IAbxlJ\nAY8fPkWTAU9ddNFFatOmjdatW6eqqqp625wmGjk5Odq7d68mTZokqXrIQ05Ojuu1tzSzZs3yugRf\niUQi6t27t2677TbdeOONXpfjGzNmzNDVV1/d6ufwfPjhh3rzzTfVu3dvXX/99Vq5cqVmzpzpdVme\nq/1descdd2jq1Klas2aNxxUB8CuaDHguNTVVffv21bp165LO0YhGo9qxY8dpNRpjxozR9u3btW/f\nPlVVVenVV1/V5MmTz0XpLcrhw4e9LsE3rLW68847FQgEdO+993pdjucOHjyoo0ePSpIOHDigZcuW\naeDAgR5X5a3HHntMRUVFKiwsVO/evXXNNddo0aJFXpflqYqKCpWVlUmSSktL9e6773LHOqB2nQwm\nfjdAkwFf6Nq1qwYMGKBQKKRQKFRvW+3tbffu3atIJNJkoxEMBvXMM89o+/btGjRokKZPn64hQ4ac\ny/J975ZbbtFnn32mbdu2qWfPnnr++ee9LslTH374oV566SWVlZVp5MiRGjlypN5++22vy/JMSUmJ\nrr76ag0fPlxbt27VxIkTuWUrGigtLdUVV1yhESNG6LPPPtPXv/71usQYAE7l4/4HrU2nTp3Url07\nbdiwoUGiUXunm1gspj179ujiiy9u9O43N9xwg4YOHcotFmu88sor2r59O+ejxhVXXCFrrfLz8zkn\nkoYPH64NGzZIkvLz87VgwQKPK/KXzMxMLVmyxOsyPNe3b19t2rRJUvXnZP78+R5XBMDPaDLgK4FA\nQMOHD9dHH32kY8eONdhem2iMGpWv8vKG2xPt3zoFJMUavNt6z8ep/vf8cE4aqj4niT9DrVXDzwnn\np6n/dzIzO+j48aPNVA3gkdrhUmiA0wLfad++vdLT0/Xpp582uOtU7V9q1Q3GwuYvrsVYKOkRr4vw\nsQcl/cTrInzufvEZasyD4ndQ48rKFnpdAgAPMScDvpSSkqK8vDxVVlbqwIEDXpcDAACQGLewTYgk\nA76Vlpam9PR07dq1K+E6GgAAAPAnkgz4mjFGeXl5Ki4ubrCOBgAAAPyJJgO+FwwGNXr06Lq1MgAA\nAHyBdTKSoslAi5CSkqL09HSFQiGFw2GvywEAAEAjaDLQorC6LAAA8A2SjKRoMtCiGGPUtm1br8sA\nAABAI3zc/wAAgJYsPz9fktSlSxctXbrU42oANCeaDAAAcE4UFBQk3da7d29lZmYqEAgoGAw2ui/g\nW6z4nRSnBQAAeOK9995Tly5dvC4DwDlAkwEAAAA45eNVt73ExG8AANDsjDG67rrrlJeXp+eee87r\ncgC4jCQDAAA0uw8++EA5OTk6cOCAJk6cqIEDB+rKK6/0uiwALiHJAAAAzS4nJ0eS1K1bN02dOlVr\n1qzxuCLAAdbJSIomAwAANKuKigqVlZXVPX/33XdZbBU4z/i4/wEAAOej0tJSTZ06VZIUjUZ16623\natKkSR5XBTjALWyT4rQAAIBm1bdvX23atMnrMgCcQwyXAgAAAOAqkgwAAADAKdbJSIgkAwAAAICr\naDIAAAAAuIrhUgAAAIAT3F0qKZIMAAAAAK6i9wIAAACcIMlIiiQDAAAAgKtoMgAAAAC4ioAHAAAA\ncMKIdTKSIMkAAAAA4CqSDAAAAMAJJn4nRZIBAAAAwFU0GQAAAABcRcADAAAAOMXVdEIkGQAAAABc\nRe8FAAAAOMHE76RIMgAAAAC4iiYDAAAAgKsIeAAAAAAnWPE7KZIMAAAAAK4iyQAAAACcYOJ3UiQZ\nAAAAAFxFkwEAAADAVQQ8AAAAgFMt+GraGNNJ0muSeksqlDTdWnskwX4XSfq1pFxJVtIN1trCxo5N\nkgEAAAC0Tg9IWmGt7S9pRc3rRF6U9FNr7SBJYyUdaOrALbj3AgAAADzU8m9hO0XSVTXPfyvpz5Lu\nP3kHY8xgSUFr7TJJstaWn86BSTIAAACA1qm7tbak5vl+Sd0T7HOJpKPGmD8aYzYYY35qjGmytSLJ\nAAAAAFquLsaYgpNeP2etfa72hTFmuaQLE3zd/JNfWGutMcYm2C8o6SuSRknao+o5HLMlPd9YUTQZ\nAAAAgBP+WCfjkLU2P9lGa+2EZNuMMaXGmGxrbYkxJluJ51oUSdpord1Z8zWLJY1XE00Gw6UAAACA\n1ulNSbNqns+S9EaCfdZK6miM6Vrz+hpJW5o6ME0GAAAA0Do9LmmiMWa7pAk1r2WMyTfG/FqSrLUx\nSf9H0gpjzCeqzm9+1dSBvQ94AAAAgJbIH8OlHLPWfinp2gTvF0iae9LrZZKGn8mxSTIAAAAAuKoF\n914AAACAx1r2OhnnDEkGAAAAAFfRZAAAAABwFcOlAAAAACda+MTvc4kkAwAAAICr6L0AAAAAJ0gy\nkiLJAAAAAOAqmgwAAAAAriLgAQAAAJxguFRSJBkAAAAAXEXvBQAAADjFit8JkWSg2Rw/flzRaNTr\nMgAAAHCO0WSg2aSmpqqqqkrr1q1TeXm51+UAAADgHGG4FJpNu3btlJ6erosvvliffvqpwuGwqqqq\nvC4LAADAGSZ+J0WSgWbXsWNHjR07VoFAQGvXrlVlZaVisZjXZQEAAMAlNBnwhDFGqampGj9+vCRp\n9erVKikp8bgqAACAM1CbZHj58CmaDHgqEAgoLS1NY8aM0ZEjR1RRUaEjR454XRYAAADOAk0GfKFN\nmzYaPHiw2rVrp507dyoUCikUCnldFgAAABzwcciC1iglJUV5eXl6//33tWnTJoXDYUUiEaWmpp6y\nZ0DSQg8qbClSJD3odRE+liLpfq+L8Dk+Q41LEb+DmsLiAWgl+KgnRJMBXwoGgxo/frz+8pe/aM2a\nNcrNzT1lj5j0kPWkthbhYSP9K+cnqflGepzz06gH+Aw1ar7hd1BTHjZeVwDAQwyXgm8ZY9SmTRuN\nGzdOlZWVKi8v18GDB70uCwAAAE0gyYDvBYNB9e/fX6WlpSopKVFFRYXXJQEAALBORiNIMtBipKSk\naPjw4Wrbtq3XpQAAAKAR9F5ocQIBZlgBQEuQn58vSerSpYuWLl3qcTXAOUCSkRSnBQAAnBMFBQVJ\nt8ViMeXn5ysnJ0dLlixpxqodhoQ6AAAgAElEQVQANAeGSwEAgGb39NNPa9CgQV6XAeAcockAAADN\nqqioSG+99Zbmzp3rdSnA2akdLuXlw6doMgAAQLO655579MQTTyglhcsQ4HzF/90AAKDZLFmyRN26\ndVNeXp7XpQCusAFvH35FkwEAAJrNhx9+qDfffFO9e/fWjBkztHLlSs2cOdPrsgC4jCYDAAA0m8ce\ne0xFRUUqLCzUq6++qmuuuUaLFi3yuiwALvPxdBEAAADAv6yRYlxNJ8RpAQAAnrjqqqt01VVXeV0G\ngHOAJgMAAABwgiQjKeZkAAAAAHAVTQYAAAAAVxHwAAAAAA5YI0UDXv+bfdzj75+Y12cFAAAAwHmG\nJAMAAABwwBqjWNDry+kqj79/YiQZAAAAAFxFkwEAAADAVV7nOwAAAECLFQsEvC7Bl0gyAAAAALiK\nJAMAAABwwMooJpKMREgyAAAAALiKJgMAAACAqxguBQAAADhgZRRluFRCJBkAAAAAXEWTAQAAAMBV\nDJcCAAAAHIpxOZ0QSQYAAAAAV9F6AQAAAA6wTkZyJBkAAAAAXEWTAQAAAMBVDJcCAAAAHGC4VHIk\nGQAAAABcRZIBAAAAOESSkRhJBgAAAABX0WQAAAAAcBXDpQAAAAAHrIyiDJdKiCQDAAAAgKtIMgAA\nAAAHqm9hy+V0IiQZAAAAAFxFkwEAAADAVeQ7AAAAgEOsk5EYSQYAAAAAV5FkAAAAAA5UT/wmyUiE\nJAMAAACAq2gyAAAAALiK4VIAAACAA1Zixe8kSDIAAAAAuIokAwAAAHCEFb+TIckAAAAA4CqaDAAA\nAACuIt8BAAAAHGCdjORIMgAAAAC4iiYDAAAAgKsYLgUAAAA41JKHSxljOkl6TVJvSYWSpltrjyTY\n7wlJX1d1QLFM0vestbaxY5NkAAAAAK3TA5JWWGv7S1pR87oeY8xlki6XNFzSUEljJH21qQOTZAAA\nAAAOnAcTv6dIuqrm+W8l/VnS/afsYyW1ldRGkpGUKqm0qQOTZAAAAAAtVxdjTMFJj789g6/tbq0t\nqXm+X1L3U3ew1q6S9J6kkprHn6y1W5s6MEkGAAAA0HIdstbmJ9tojFku6cIEm+af/MJaa40xDeZZ\nGGP6SRokqWfNW8uMMV+x1v61saJoMgAAAAAHrIyiPh8uZa2dkGybMabUGJNtrS0xxmRLOpBgt6mS\nVltry2u+5h1Jl0pqtMlguBQAAADQOr0paVbN81mS3kiwzx5JXzXGBI0xqaqe9M1wKQAAAOBcibXs\ny+nHJf3OGHOnpN2SpkuSMSZf0nestXMl/UHSNZI+UfUk8KXW2v9u6sAt+qygZamqqlI8Hpe1VsYY\nr8sBAABo1ay1X0q6NsH7BZLm1jyPSbr7TI9Nk4FmEw6HFQ6H9dFHHykYDCocDquoqEixWEzxeFwp\nKYzeAwAAOB/QZKDZZGVlKT09XZdddpkikYhWrVqlaDSqqqoqffzxx7LW6sSJE9q1a1fd+wAAAH51\nHqyTcc7QZMATqampCgaD6t27t4qLi3XppZcqHo/rww8/VJs2bRSNRrVx40aVl5dr3bp1ysrKUiQS\nUSgUqj5AIFV6mCFXSaUEpfmcn6RSgtIDnJ9G8RlqXEqQ30FNCaR6XQEAD9FkwDdSUlIUCASUk5Oj\n3bt3a+zYsfroo480ePBglZWVqaioSNu2bVN5ebkUi+g++7DXJfvWk+YhLbSnLtiJWgvNT/Swvc/r\nMnztIfMkn6FGLDQ/4XdQE540D3ldAnDOkWQkxyB4+F67du3UrVs3paWladSoUWrfvr3XJQEATkN+\nfr7y8/M1adIkr0sB0MxIMgAAwDlRUFCQ8P1wOKwrr7xSlZWVikajmjZtmn70ox81c3UAziWaDAAA\n0KzS0tK0cuVKtW/fXpFIRFdccYW+9rWvafz48V6XBpwxv6/47RWGSwEAgGZljKkb+hqJRBSJRFg/\nCTjPkGQAAIBmF4vFlJeXpy+++ELz5s3TuHHjvC4JOGPVE7+5nE6EJAMAADS7QCCgjRs3qqioSGvW\nrNHmzZu9LgmAi2gyAACAZzp27Kirr75aS5cu9boUAC6iyQAAAM3q4MGDOnr0qCTpxIkTWrZsmQYO\nHOhxVcCZq10nw8uHXzGIDAAANKuSkhLNmjVLsVhM8Xhc06dP1ze+8Q2vywLgIpoMAADQrIYPH64N\nGzZ4XQaAc4gmAwAAAHDIz0OWvMScDAAAAACuIskAAAAAHLAyrPidBEkGAAAAAFfRZAAAAABwFcOl\nAAAAAAeq18ngcjoRkgwAAAAArqL1AgAAABziFraJkWQAAAAAcBVNBgAAAABXMVwKAAAAcKB64jfD\npRIhyQAAAADgKpIMAAAAwAGSjORIMgAAAAC4iiYDAAAAgKsYLgUAAAA4FGW4VEIkGQAAAABcRZIB\nAAAAOFA98ZvL6URIMgAAAAC4iiYDAAAAgKvIdwAAAAAHWCcjOZIMAADgewsWLNBTTz1V93r+/Pl6\n+umnPawIQGNoMgAAgO/NmTNHL774oiQpHo/r1Vdf1cyZMz2uCpBiCnj68CuGSwEAAN/r3bu3Onfu\nrA0bNqi0tFSjRo1S586dvS4LQBI0GQAAoEWYO3eufvOb32j//v2aM2eO1+UAaARNBgAAaBGmTp2q\nBQsWKBKJ6OWXX/a6HEBWhhW/k6DJAAAALUKbNm109dVXq2PHjgoEuLAD/IwmAwAAtAjxeFyrV6/W\n73//e69LAdAEmgwAAOB7W7Zs0Te+8Q1NnTpV/fv397ocQFLtOhlcTifCWQEAAL43ePBg7dy50+sy\nAJwmmgwAAADAIT+vVeElFuMDAAAA4CqaDAAAAACuYrgUAAAA4ED1xG+GSyVCkgEAAADAVSQZAAAA\ngEMkGYmRZAAAAABwFU0GAAAAAFcxXAoAAABwwMooynCphEgyAAAAALiKJAMAAABwoPoWtlxOJ0KS\nAQAAAMBVNBkAAAAAXEW+AwAAADjEOhmJkWQAAAAAcBVJBgAAAOBA9cRvkoxESDIAAAAAuIomAwAA\nAICrGC4FAAAAOMCK38mRZAAAAABwFUkGAAAA4BArfidGkgEAAADAVTQZAAAAAFxFvgMAAAA4wDoZ\nyZFkAAAAAHAVTQYAAAAAVzFcCi1KNBpVLBbzugwAAACGSzWCJAO+FQ6HFY1GtWPHDm3cuFHl5eUq\nKChQVVWV16UBAE5Dfn6+8vPzNWnSpHrv7927V1dffbUGDx6sIUOG6Omnn/aoQgDnCkkGfKG8vFzH\njx9XOBzW2rVrVV5ers2bNysajSo9PV3du3dXKBTS+PHj9dFHH3ldLgDgNBQUFCR8PxgM6sknn9To\n0aNVVlamvLw8TZw4UYMHD27mCoGzR5KRGE0GmpW1VkeOHFFZWZlOnDihVatWqby8XDt27FBmZqaC\nwaCGDx+udevWKT8/Xx999JGys7MbHCclNUVPmoc8+AlahpRgihaan3hdhm+lBFP0kHnS6zJ8jc9Q\n41KC/A5qSkpq8sES2dnZdb/bMzMzNWjQIO3bt48mAziP0GSg2Xz55ZcKhUIqLi5WVlaW2rRpo7Fj\nx+rjjz/WiBEjJEn79+9XWlpak8eKR+Kyr5zrilsuc0tcltEHSZnvxWXpMRpl7uMz1BjzPX4HNcXc\nEj+t/QoLC7VhwwaNGzfuHFcEoDnRZKDZdO7cWRkZGRoyZIik6jG5gQARIwC0VuXl5brpppv01FNP\nKSsry+tygDNmZRRluFRCTPwGAADNLhKJ6KabbtJtt92mG2+80etyALiMJgMAADQra63uvPNODRo0\nSPfee6/X5QCOVd/CNujp42wYY75tjPnUGBM3xuQ3st8kY8w2Y8wXxpgHTufYNBkAAKBZffjhh3rp\npZe0cuVKjRw5UiNHjtTbb7/tdVlAa7RZ0o2S3k+2gzEmIOlZSV+TNFjSLcaYJu/SwJwMAADQrK64\n4gpZa70uA2j1rLVbJckY09huYyV9Ya3dWbPvq5KmSNrS2BfRZAAAAAAO+WCdjC7GmJMXpXnOWvuc\ni8fPkbT3pNdFkpq8HRxNBgAAANByHbLWNjafYrmkCxNsmm+tfeNcFUWTAQAAADhQPfHb8ySjUdba\nCWd5iH2Sck963bPmvUYx8RsAAABAMmsl9TfG9DHGtJE0Q9KbTX0RTQYAAADQChljphpjiiRdKukt\nY8yfat7vYYx5W5KstVFJfy/pT5K2SvqdtfbTpo7NcCkAAADAgZa+4re19nVJryd4v1jSDSe9flvS\nGd1nmiQDAAAAgKtIMgAAAACHznbV7fMVSQYAAAAAV9FkAAAAAHAV+Q4AAADgQEtYJ8MrJBkAAAAA\nXEWTAQAAAMBVDJcCAAAAHGC4VHIkGQAAAABcRZIBAAAAONSSV/w+l0gyAAAAALiKJgMAAACAqxgu\nBQAAADhQPfGby+lESDIAAAAAuIrWCwAAAHCAW9gmR5IBAAAAwFU0GQAAAABcxXApAAAAwCGGSyVG\nkgEAAADAVSQZAAAAgANM/E6OJAMAAACAq2gyAAAAALiK4VIAAACAA1ZSlOFSCZFkAAAAAHAVSQYA\nAADgiFGMy+mESDIAAAAAuIomAwAAAICryHcAAAAAB1gnIzmSDAAAAACuIskAAAAAHCLJSIwkAwAA\nAICraDIAAAAAuIrhUgAAAIADVoYVv5MgyQAAAADgKpoMAAAAAK5iuBQAAADgQPU6GVxOJ0KSAQAA\nAMBVtF4AAACAQ6yTkRhJBgAAAABX0WQAAAAAcBXDpQAAAAAHqid+M1wqEZIMAAAAAK4iyQAAAAAc\nsDKKxUkyEiHJAAAAAOAqmgwAAAAArmK4FAAAAOCElaJRhkslQpIBAAAAwFUkGQAAAIAD1hrFolxO\nJ0KSAQAAAMBVNBkAAAAAXEW+AwAAADhQPVyKid+JkGQAAAAAcBVJBlqUSCSicDjsdRkAgNOQn58v\nSerSpYuWLl1a9/6cOXO0ZMkSdevWTZs3b/aqPODsWZFkJEGSgRYhHo+rqqpKa9asUUoKH1sAaAkK\nCgpUUFBQr8GQpNmzZzd4D8D5has1+N6BAwe0evVqWWs1btw4tWnTxuuSAABn4corr1SnTp28LgPA\nOcRwKfjW8ePHVVFRodLSUo0ePVrr169XMMhHFgAA+IO1RtEIw6US4YoNvhMOh3XixAlt27ZNbdu2\n1bBhw7wuCQAAAGeAJgO+EY1GVVlZqfXr1ys1NVX5+flatWpVwn1TA5K5pZkLbEGCKZL5ntdV+Fcw\nRTL3eV2Fv/EZalyQ30FNSuUfd9EqGMVjXE4nwlmBLxQVFWn37t0yxmj8+PFavXq1jDFJ94/EpF81\nY30tzV1xaaHXRfjYQs5PkzhHjVvI76Am3RXzugIAXmLiNzz15ZdfqqKiQuXl5Ro7dqzatGnD3aMA\n4Dx3yy236NJLL9W2bdvUs2dPPf/8816XBMBlJBnwREVFhUKhkPbs2aN27dpp4MCBXpcEAGgmr7zy\nitclAO6wklgnIyH+yRjNylqrrVu36pNPPlFaWppGjRpFcgEAAHCe4eoOzebw4cOqqKhQhw4dNG7c\nOAUCdP4AAADnI4ZLodlkZWUpIyNDPXr08LoUAACAs2cNw6WSIMlAswkGg43eMQoAAADnB5IMAAAA\nwAkrKco/oCZCkgEAAADAVTQZAAAAAFzFcCkAAADAqajXBfgTSQYAAAAAV5FkAAAAAE5YkWQkQZIB\nAAAAwFU0GQAAAABcxXApAAAAwAmGSyVFkgEAAADAVSQZAAAAgBNWUsTrIvyJJAMAAACAq2gyAAAA\ngFbIGPNtY8ynxpi4MSY/yT65xpj3jDFbavb93ukcm+FSAAAAgBNWUszrIs7KZkk3Svr/GtknKuk+\na+16Y0ympHXGmGXW2i2NHZgmAwAAAGiFrLVbJckY09g+JZJKap6XGWO2SsqRRJMBAAAAnBOt6Ba2\nxpjekkZJ+ripfWkyAAAAgJarizGm4KTXz1lrn6t9YYxZLunCBF8331r7xul+E2NMe0n/T9I91trj\nTe1PkwEAAAC0XIestQknbUuStXbC2X4DY0yqqhuM/2ut/ePpfA1NBgAAAOBEK1jx21RP2Hhe0lZr\n7c9O9+u4hS0AAADQChljphpjiiRdKuktY8yfat7vYYx5u2a3yyXdLukaY8zGmscNTR2bJAMAAABw\nooUnGdba1yW9nuD9Ykk31Dz/QFLy208lQZIBAAAAwFU0GQAAAABcxXApAAAAwIkWPlzqXCLJAAAA\nAOAqmgwAAAAArmK4FAAAAOAEw6WSIskAAAAA4CqSDAAAAMApkoyESDIAAAAAuIomAwAAAICrGC4F\nAAAAOGElRbwuwp9IMgAAAAC4iiQDAAAAcMJKinldhD+RZAAAAABwFU0GAAAAAFcxXAoAAABwghW/\nkyLJAAAAAOAqkgwAAADACZKMpEgyAAAAALiKJgMAAACAqxguBQAAADjBcKmkSDIAAAAAuIokAwAA\nAHCKJCMhkgwAAAAArqLJAAAAAOAqhksBAAAATjDxOymSDAAAAACuoskAAAAA4CqGSwEAAABOMFwq\nKZIMAAAAAK4iyQAAAACcsJIiXhfhTyQZAAAAAFxFkwEAAADAVTQZaHGstV6XAAA4Dfn5+crPz9ek\nSZMabFu6dKkGDBigfv366fHHH/egOsAFVlLM44dPMScDLUo4HFYoFPK6DADAaSgoKEj4fiwW07x5\n87Rs2TL17NlTY8aM0eTJkzV48OBmrhDAuUKTgRYjHo9r/fr1atu2rdelAADOwpo1a9SvXz/17dtX\nkjRjxgy98cYbNBlombiFbUIMl0KLUFZWplAopKFDhyoQCHhdDgDgLOzbt0+5ubl1r3v27Kl9+/Z5\nWBEAt5FkwPeOHj2qTz/9VOnp6crKyvK6HAAAADSBJAO+Fo1GtWXLFo0ePVopKXxcAeB8kJOTo717\n99a9LioqUk5OjocVAQ7Vrvjt5cOnuGqDb5WWlqqyslJ5eXlq166d1+UAAFwyZswYbd++Xbt27VJV\nVZVeffVVTZ482euyALiI4VLwpUgkosLCQqWnpystLc3rcgAADuTn59c979Kli5YuXSpJCgaDeuaZ\nZ3T99dcrFotpzpw5GjJkiFdlAs7VJhlogCYDvrN7925FIhFdccUV+vjjj70uBwDgULJb2ErSDTfc\noBtuuKEZqwHQnGgy4BvWWlVWVurIkSNKT09PeBepeDyulJQUtW/XTnedOOFBlS1DiqSFXhfhY5yf\npnGOGheQdJfXRfhcQNKwYcMapNEnJxoAzl80GfAFa60+++wzxeNxjRgxQqtWrUq4XzQa1ejRo1Vc\nWlrva7dt26ZgMKiLL75YxpjmKtv3Dh8+rOLiYg0dOtTrUnwlHo+roKBAY8eO9boU31m/fr2GDRum\n1NRUr0vxlc8//1xZWVm68MILvS7FV/bs2aNjx45pyJAhDW7OsW7dOv3DP/yDXn75ZV1yySUeVQic\nY1ZSxOsi/IkmA76wefNmpaamql27dgmbhNoEY+3atQ22hcNhSVLbtm0bjeZbG2utQqGQ2rVrl/C8\ntWaxWExVVVWclwTC4bAKCgoUDPLXw8mstSouLtaePXv4h4xTVFVV6YMPPlDbtm3rnZsHHnhA4XBY\nl112mbp3797gBh4kGsD5jb9F4Kl4PK5QKKTu3burb9++CROMeDyuSCSi0aNH13vfWqvt27erY8eO\n6tevH3/xn2LPnj2KxWLq06eP16X4TklJiSKRiC666CKvS/Gd/fv3q7KyUr169fK6FN8pLi5WKBRS\nv379vC7Fd3bv3q3y8nINHjy47nfxihUrJFX/I9Lf/u3f6te//rVGjhzpZZmA+6ykmNdF+BO3sIVn\nrLVav359k8OcotFog23WWn3xxReSRIORQGVlpfbv389FdBJlZWXKzMz0ugxfyszMVFlZmddl+FJ2\ndraOHTumiooKr0vxnV69eikjI0Nbt26VtbbetqFDh+qFF17Q7bffTtoMtCIkGfBEJBJRKBRSnz59\nVFhYmHCfxoZIVVZWylrLEKkkTpw4odTUVK1fv97rUnwpFArp6NGjNKdJVFRUMJQsiVgspoKCAqWn\np3tdii9VVlbqr3/9a72hUQ888ICOHTumeDyuSZMmqWvXrsrIyKj3dQydAs4/NBlodpWVlVq/fr3S\n0tLUo0ePhE2GtbZukvep7+/cuVORSEQDBgzgIjGBY8eOqbCwUMOHD+f8JGCt1dq1a5n03Yh169Zp\n+PDhTP5OYuvWrerSpYu6du3qdSm+Y63Vrl27VFlZqYEDB8oYUzdsSpJ27typ22+/XU8//bQuv/xy\nDysFXMQ6GQnRZKBZxeNxrVu3TgMGDNC2bdsS7hOLxRQIBEgwHKqoqFC7du04P0nE43FVVlbyL/WN\nYPJ346y1Ki0tVXp6Oo18EpWVlXWTwWvVJhrWWk2ZMkWdO3duMGyRRAM4f/A3CJpNeXm5QqGQxo0b\npw4dOjTYbq2VtVaxWKxBgiFJu3bt0okTJzRo0CD+Yk9i3759CofDuvjii70uxbeY2Nw0JsY3be/e\nvYpGo9xYIQlrrXbs2KFYLKZLLrmkQaJRVFSkW2+9VY899pgmTpzoYaXAWWLF76RoMtBsYrGY2rVr\nl7TBiMfjjSYY8Xicf6FvRO0ta9PT03X48GGvy/GtcDisYDCoAwcOeF2Kb9Xe4rf0pPVoUF/t/28H\nDhxosD4E/lc4HNaBAweSJhq33nqrOnbsmPDvBVINoGWjyUCz6dChQ8JVvKXqISzJEoxEt0ZEQ9u2\nbVOvXr3UvXt3r0vxNRaba1rtYoVjxozxuhRfO3LkiIqKijRs2DCvS/Gt2luNG2Pq7gR4cqJx8OBB\nTZ8+XQ899JCmTJniYaUA3MY/v8BTJw+RSmTPnj0qKytjiFQTysrKFAqF1K1bN69L8bXaGwrQYDQu\nJSVFgUBA0ShjABpzwQUXyBijL7/80utSPHX48GFNmTJFI0eO1JQpU3TkyJG6bcYY9e/fX/F4XP9/\ne/ceXGV17nH8t7JzI4FcQCAJSQghgBKsQBOQOoSKQQI5RpFCSQtjRUZEkFpGMBQKtYVDLFPkPkUR\nhFMUhApMAw2CIx28RUsVJIFDkIC5kasRczO39/yBydFKuGxf2Yl8PzP5g2TvN4s1w5Anv+dZa+XK\nlRo4cKAGDhyobdu2SZK6du2qxYsXKzk5WUFBQZo9e3bLEbg7d+5suUn8PxPsZcuWKSoqSv369dOB\nAwdaPp+enq5+/fopKipKqampN+Bvj5te843frvxoo0gy4FJXapGqq6tTY2OjvL29dfToUResrv2o\nrq6Wl5cXrWRXwdD3tWse/m4tfcQlTU1N+vjjj2/qIfANGzaod+/eWrhwoV555RXNmzdPjz766Dde\nc/HiRa1fv16rV6+Wl5eXZsyYoeDgYC1dulTHjx9Xt27dVFNTo40bN+of//iH/P39VVtbK2OMAgIC\nvvGsrKwsbd++XZmZmSooKFB8fLxOnz4tSZo5c6YOHjyo0NBQxcbGKikpSf37979hewHg/1FkwCWu\nNuSdl5en8vJyDRgwgH7nq7hw4YIuXryovn37unopbV5RUZFqamoUERHh6qW0eQUFBWpsbFRYWJir\nl9LmnTt3Tm5ubjftoPz06dO1f/9+eXp6KisrSx988IGKi4v10ksvKTAwUNKlVGLMmDGKiIhQenq6\nKisr9eijj2rOnDmqqKjQ0aNHdeTIEY0bN055eXkaM2aMVq1aJWOM7rrrLs2cOVPl5eWKiIjQsGHD\nNGnSJK1evVrbtm1TSUmJ+vTpo/Pnz2vEiBGKjIxURESEamtrFR8fr5CQEH4BA7gARQZuuObiggTj\nu/v6sDe/nb+6L7/8Ug6HQyUlJa5eSpvXPPx94cIFVy+lzWv+d1hYWHhT/lKksLBQubm52rBhg267\n7TYdO3ZMvXv3/kai8f7778sYo4KCAq1fv14JCQny9fXV8uXL1aNHD91zzz06fvy4AgICVFFRoRde\neKEl0Wj+fHh4uM6fP6+ysjLNmTNHkydP1ty5c/XII4+oa9euSktLU2RkZMu6Fi1apKysLK1du9ZV\nW4ObgSXp8h3fNz2KDNxwjY2NampqumyCkZ+fr9LSUt1+++035X/W1+vMmTMKDQ1VSEiIq5fSLnz4\n4YeKjo6Wp6enq5fS5jXfacPw97UpLS1VUVGRoqOjXb2U70VSUtJlTxtbtGiR3N3dFRsb25Jo7Ny5\n81uJRlhYmGpra1VeXq4xY8bI09NT58+fV21trU6ePKk//vGPqqio0Jo1a7Ro0SIdO3ZMeXl5SkhI\nUGZmprZu3aohQ4bogQce0Ntvv63U1FQlJiYqMDBQhYWF2rJliwICArR79+5vJErnz59Xv3791NjY\nqGnTpiklJeVGbhtwU6PIwA3T1NSkpqYm1ddffkqpuLhYJSUlio6OvuIwOC6pqalRRUWFIiIi2Ktr\nYFmW6urq5HA42K/r0LxnuLLAwEDl5eXps88+k5+fn6uXY7vdu3e3+rWuXbsqPz9fxcXFamxslJub\nm+Lj43XixAkNHz5cf/7zn3X77bfrrbfekqenp3r06KFPPvlEGRkZuv/++5Wenq5Vq1apV69eys3N\n1SeffKJnnnlGy5cv16FDh1RXV6dbbrlFqampGjNmjD744AOFhoYqNTVVzz77bEvxs23bNq1YsUKL\nFi3S5s2btWTJEhUXF+uZZ57R3LlzmdHA94N7Mlplmk9xuEbX9WLg6+rr6/Xuu++2+gNeU1OTjDE3\n7fCkMyzLYr+uA/t1fdiv69P8/+kPec+eeuqpb93DU15erqFDh+qdd95RcnKyduzYoU2bNunhhx/W\n5s2b9Zvf/EZr167VY489poSEBNXX1ystLU1Dhw5VSEiI9u3bp+DgYF28eFEBAQHKzMxUUFCQysvL\n5ePjo5KSEt122206e/as+vXrpxMnTqhDhw4KCQnRgQMHNGjQIHl5eSk/P199+/ZtOSJ31KhR6ty5\ns7788kutWbNGb7/9tmOd+xUAABOOSURBVCRp/vz5N3zf4JR28Q/JdImxlOjimZ//MUcty4px7SK+\njSIDAAA4raysTBMnTtSRI0c0dOhQHT9+XKdOndKdd96p4cOHa8eOHYqIiJCnp6eKi4tVV1enhIQE\n3X777crNzVV4eLiWLVumqqoqdezYUeHh4Vq4cKE2b96sDz74QBUVFXJzc5NlWcrLy9Pdd9+t8ePH\na9myZfLw8Ghpf+zdu7cuXryo/Px8SZKfn5+GDx+uO+64Qx07dtSzzz6rpqYmhYWFyd3dnWHwto8i\n41q10SKDpncAAHBd4uPjNWDAAA0YMEAjRoxQUVGR/Pz81KtXLxljtGXLFv385z9XSEiI3N3dlZ2d\nrSlTpig5OVmdO3dWdHS0ampq9Prrr2vSpEny9PRUSEiIli1bprNnz+rFF1/Uk08+qaCgIIWEhGjJ\nkiVyOBxasGCB7r//fjU0NMjb21slJSVqbGzUr371Kx07dkwzZ85Ujx49VFZWpunTp+utt97S66+/\nrgEDBkiSHnjgAX300UcUGLBXg4s/2ihmMgAAwHU5dOjQtz7XnGhUV1dr37592rt3r6Kjo9XQ0CA/\nPz/NmjVLBw4c0O9+9zvNnz9ftbW1Wr58uSZPnqzKykr5+vpqxYoVqqurk5ubm44cOaLs7Gz5+vpq\nx44dsixL6enpWr16tX75y1/KsixNmzZN9957rw4ePKicnBzt3LlTubm5Cg0NVWBgoMrKyjRr1ixt\n3rxZ5eXlSk9P17lz5zjGGrgBSDIAAMB31qVLF73xxhv69a9/rcTERPn7+6uoqEgPP/ywSktLtW/f\nPhUUFGjq1Kk6efKkOnfurIsXL6qhoUG+vr5avHixYmNjFRwcrK5du2rXrl3y9/fXc889py+++ELj\nx4+Xn5+fUlJSZIzR7Nmzdfr0ac2ZM0dFRUV6+umndd999ykkJESTJk3SoEGD5OvrK2OM/Pz8Wg4w\nGDRokJ5//nkX7xZ+MLjxu1XMZAAAAKfFx8d/4z6VxsZG5ebmysvLS5999plKS0v1/vvvKzk5WRUV\nFVq2bJlSUlK0adMmzZo1Sx4eHvrpT38qLy8v/e1vf1OXLl1UVVUld3d3+fn5acaMGVqyZIk6duyo\nqqoqWZal+vp6vfzyy3r88cdVXV2t6upqORwO/fjHP5a/v7+OHTumkpISxcXF6f3331dgYKDi4uK0\nf/9+VVVVycPDQ8nJydq0aZMLdw5X0T5mMjrHWLrHxe13u5jJAAAAPzCHDh3SiRMnWj5OnjypyspK\nbdiwQZ06dVJ1dbVmzpyp2bNny8/PT6+88oqysrI0depURUZGaunSpSorK1Pnzp0VERGhW2+9VZ06\nddLSpUvl7u6uHTt2aOLEiUpMTJRlWZo0aZJ69eqluXPnKiUlRePHj2+Z6Wg+FtfLy0u9e/dWcHCw\nAgIClJGRoenTp2v48OHq2bOnnnjiCY6yBb5nFBkAAOB70bNnTy1dulRRUVE6efKkevfurUmTJmnv\n3r0trxkxYoSOHj2qcePGqampSTk5OaqoqJC/v7/q6+v1k5/8RG5ubvroo48UFBSkkpISxcbGytPT\nU926dVN2drYGDx6sqqoqffrpp3ruuee0atUqeXl5aciQIaqvr29JO6RLRw2/+eabLcPgwHfSfOO3\nKz/aKNqlAACA7d59910tWLBAJSUl+vTTT9W1a1clJyertrZWf/nLX2SMUVBQkLZs2aKf/exncnd3\nV15entavX6/Zs2e3zFN4eHjo4sWLGjt2rDIzM3X27FkFBga2XOwaGBiohIQErV+/Xh07dlR9fb1C\nQ0NVWFio2tpaORwOBQcHq7CwUA0NDWpqalJUVJSGDBmijIwM+fj46KWXXtLgwYNdvGP4D+2jXSow\nxtLdLm6X2k27FAAAuEnExsbq3LlzeuyxxzR+/Hj5+vpqwoQJ2rp1qx544AGVlpaqurpaK1euVGBg\noObNm6fx48ersbFRHh4eCg4O1uDBgxUcHKygoCAtXrxYRUVF8vf314wZM1RfX68ePXrotdde086d\nOxUWFqb169erb9++Gjx4sO644w717dtXU6ZMUUJCgkpKSvTiiy9qwoQJmjx5svbs2aPs7Gw9//zz\nmjFjhqu3C+1V843fHGH7LRQZAADAdu7u7lq7dq1SU1O1fft2TZw4UVVVVaqvr5cxRp6enpo+fbqy\nsrJ09uxZrVu3TqmpqdqzZ48cDoeSk5P1r3/9S8eOHdPatWuVlpamLl266Omnn9bGjRtVVVWlYcOG\nqaqqSv7+/po4caJWr16tvLw8HT58WF26dFFsbKxGjhypsrIyDR48WOvXr1dqaqoKCgrk4+OjsrIy\n3XnnnaqoqFBhYaGrtwz4QaFdCgAAfG8aGhrUt29fvfHGG8rIyNDjjz+uI0eOKDo6WnPnztXGjRvl\n7u6uyMhIZWRkKCwsTHfccYfS0tL04YcfKiYmRmFhYaqsrJRlWSouLpZlWfL19ZWfn588PDz0xRdf\n6L333lN0dLRGjx6tN954Qz179lRAQIAqKyvl5eWl4uJide/eXcYYZWVlyRij2tpa/fOf/1R8fLwi\nIyPl4+OjBx98UIsWLXL1tqG9tEsFxFiKc3G71N/bZrsUl/EBAIDvTXOiMXr0aH3++eeKjIxUdHS0\nFi5cqBdeeEH33Xef1q1bp/DwcIWHh6usrEzz58+XdGmuw9390o8q1dXVCg8Pl8PhUFZWljw9PeVw\nOFRXV6e6ujrdeuutamxs1Ouvv66OHTvK19dXmZmZeuSRR7RixQrFxcUpJydHwcHBkqRJkybJmEs/\nxwYEBOjll19WTEyb+zkNbV1zuxS+hXYpAADwvRo7dqxOnz6tPXv2qEuXLpKkxMREdevWTf3799c7\n77wjDw8PVVVVKTQ0VG5ul3482bNnj7y9vXXq1Cn17t1bp0+f1tChQ/XSSy8pKChIu3fv1oQJE9TQ\n0KDIyEitXr1affr00W9/+1vFxMTI3d1du3btkpeXl5599lkNGTJEe/fuVUNDg9LT0zVq1Ch98cUX\nqq2tVUhIiGbPnq2oqCj96Ec/0r///W9XbhnQ7pFkAACAGyI2NlbZ2dnKycnR+fPnVVpaqsTERI0b\nN04pKSk6c+aM/v73v2vlypUaNmyYTp48qeHDh2vTpk267bbbdPr0aT300EP605/+pOrqanXs2FFv\nv/22OnTooKlTpyo1NVV1dXVKSkrSpk2b1KFDB5WVlSkiIkINDQ2aPHmyZs6cqbCwMPXv319xcXFa\ns2aNamtrddddd6myslJvvvmmKisrNWPGDGVkZLh6y9DWNd/4jW9hJgMAANww+/fv15NPPqnPP/9c\nYWFhWrNmjSZPnqykpCTV19ere/fuevXVV1VTU6OCggIdPHhQf/jDH/T73/9eiYmJ6tSpk3JzczVn\nzpyWVOTNN9/Uxx9/rPLycvXq1UunTp2SJEVFRencuXPy8vJSbGys0tLS1L17d61cuVIfffSR0tLS\ndOHCBR08eFDbtm2Tv7+/du/erezsbPXr10+HDx9uaa/CDdc+ZjL8YiwNdfFMxqG2OZNBuxQAALhh\n/rN1Kj8/X3fffbe6deumHj16KCIiQiNGjNCZM2c0fPhwubm5KT8/X8HBwXJzc9O5c+cUGBgob29v\n5efnKywsTMYYvfbaa7r33nvl5eXV8r0efPBBeXp66qmnntLhw4dVW1urmpoaTZkyRevWrdO5c+fU\noUMHxcXFKT8/v6XQKS0tVWhoqPLz8124U0D7RpEBAABuuObWqaKiIjU1NWn79u1KSkr6xmuSkpK0\nZcsWSdK+ffs0cuRIGWPk4+Oj3bt3q7GxUSUlJcrOztaQIUN0yy236MKFC8rJyVFdXZ3++te/yhij\nxYsXS5J27dold3d3eXt7S1LL4HdzV0dWVpaamppaEhLgmnDj92UxkwEAAG645lOnpk+frrKyMi1Y\nsEDR0dF64oknFBoaKkl65JFHNGXKlJZ7NNLS0tTQ0KCamhpNnz5dK1askMPh0NatW+VwOFRQUKB5\n8+Zp9OjRqq6uVl1dncLDw1sGybdv365u3bqpsLCw5RZwLy8vDRgwQBcuXNCpU6e0fft2GWOUl5en\nHj16uHKLgHaNJAMAALjE2LFjlZOTo6CgIP3iF79QXV2dysrK9PTTT0uSvL29tXPnTi1fvlxxcXGK\njIzUrl27NHLkSC1cuFBHjhxRSEiIRo4cqZycHGVnZ+vJJ5/U6dOnNWrUqJb2qoEDB2ratGk6fPiw\nkpOTW9KRLVu26KGHHlJmZqa2bt2qPn36aNiwYXrvvffk7+/PPAaurp3f+G2MmWCMyTTGNBljrjjX\nYYxxGGM+NMakXcuzSTIAAIDLfP0ejcbGRk2dOlXR0dFatGiRYmJilJSU1JJoREVFqXPnztq+fbsk\nKTo6WhMnTlT//v3l7u6udevWyeFwSJJqa2vlcDjU0NCg0tJSDRs2TJKUkpKiiRMn6sUXX1TPnj31\n6quvSrpU8Ozfv19RUVHy8fHR5s2bXbMhwI11QtKDkjZcw2t/LemkJL9reTCnSwEAAKCtaR+nS3WK\nsTTIxadLHfnup0sZYw5LesqyrMv+ZYwxoZK2SFoqaY5lWf91tWeSZAAAAADOaBs3ft9ijPl6cfC8\nZVnP2/w9VkqaJ6nTtb6BIgMAAABov0qvlGQYYw5JCrrMlxZYlrX3ag83xvyXpGLLso4aY356rYui\nyAAAAAB+oCzLiv+Oj7hLUpIxZqwkb0l+xpi/WpY1+UpvosgAAAAAnGFJqnf1Ir5flmXNlzRfkr5K\nMp66WoEhcYQtAAAAcFMyxowzxuRJGiZpnzHmwFefDzHG7P9Oz+Z0KQAAALQx7eN0KZ8YS7e6+HSp\nD7/76VLfB5IMAAAAALaiyAAAAABgKwa/AQAAAGe5/p6MNokkAwAAAICtSDIAAAAAZ7SNG7/bJJIM\nAAAAALaiyAAAAABgK9qlAAAAAGfcBDd+O4skAwAAAICtSDIAAAAAZ1iSGl29iLaJJAMAAACArSgy\nAAAAANiKdikAAADAGdyT0SqSDAAAAAC2IskAAAAAnEWScVkkGQAAAABsRZEBAAAAwFa0SwEAAADO\n4MbvVpFkAAAAALAVRQYAAAAAW9EuBQAAADjDktTo6kW0TSQZAAAAAGxFkgEAAAA4gxu/W0WSAQAA\nAMBWFBkAAAAAbEW7FAAAAOAM2qVaRZIBAAAAwFYkGQAAAIAzuPG7VSQZAAAAAGxFkQEAAADAVrRL\nAQAAAM7ixu/LIskAAAAAYCuSDAAAAMBZlqsX0DaRZAAAAACwFUUGAAAAAFtRZAAAAACwFUUGAAAA\nAFtRZAAAAACwFUUGAAAAAFtRZAAAAACwFUUGAAAAAFtRZAAAAACwFTd+AwAAAE6xJNW7ehFtEkkG\nAAAAAFtRZAAAAACwFe1SAAAAgFMsSQ2uXkSbRJIBAAAAwFYUGQAAAABsRbsUAAAA4BROl2oNSQYA\nAAAAW5FkAAAAAE5h8Ls1JBkAAAAAbEWRAQAAAMBWtEsBAAAATmHwuzUkGQAAAABsRZIBAAAAOIUk\nozUkGQAAAABsRZEBAAAAwFa0SwEAAABO456MyyHJAAAAAGArkgwAAADAKQx+t4YkAwAAAICtKDIA\nAAAA2Ip2KQAAAMAplhj8vjySDAAAAAC2IskAAAAAnMLgd2tIMgAAAADYiiIDAAAAgK1olwIAAACc\nwuB3a0gyAAAAANiKJAMAAABwCoPfrSHJAAAAAGArigwAAAAAtqJdCgAAAHAKg9+tIckAAAAAYCuK\nDAAAAAC2ol0KAAAAcAqnS7WGJAMAAAC4CRljJhhjMo0xTcaYmCu8LsAYs8sYc8oYc9IYM+xqzybJ\nAAAAAJzS7ge/T0h6UNKGq7xulaR0y7J+ZozxlORztQdTZAAAAAA3IcuyTkqSMabV1xhj/CXFSfrV\nV++pk1R3tWfTLgUAAAC0X7cYY/71tY9HbX5+L0klkjYbYz40xmw0xvhe7U0kGQAAAIBT2sTgd6ll\nWVeapzgkKegyX1pgWdbea3i+u6TBkp6wLCvDGLNKUoqk313tTQAAAAB+gCzLiv+Oj8iTlGdZVsZX\nf96lS0XGFVFkAAAAAE5r14PfV2VZ1gVjTK4xpp9lWf8r6R5JWVd7HzMZAAAAwE3IGDPOGJMnaZik\nfcaYA199PsQYs/9rL31C0jZjzHFJAyX991WfbVnW9azlul4MAAAAOKH1447aEGP6WNIKF68i6eiV\nZjJchXYpAAAAwCltYvC7TaJdCgAAAICtSDIAAAAAp5BktIYkAwAAAICtKDIAAAAA2Ip2KQAAAMAp\nln7o92Q4iyQDAAAAgK1IMgAAAACnMPjdGpIMAAAAALaiyAAAAABgK9qlAAAAAKcw+N0akgwAAAAA\ntqLIAAAAAGAr2qUAAAAAp3C6VGtIMgAAAADYiiQDAAAAcAqD360hyQAAAABgK4oMAAAAALaiXQoA\nAABwCoPfrSHJAAAAAGArkgwAAADAKQx+t4YkAwAAAICtKDIAAAAA2Ip2KQAAAMApDH63hiQDAAAA\ngK1IMgAAAACnMPjdGpIMAAAAALYylmW5eg0AAABAu2OMSZd0i4uXUWpZVoKL1/AtFBkAAAAAbEW7\nFAAAAABbUWQAAAAAsBVFBgAAAABbUWQAAAAAsBVFBgAAAABbUWQAAAAAsBVFBgAAAABbUWQAAAAA\nsBVFBgAAAABb/R8qrIiwAkp45QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3064dba860>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pp.plot_grid(g, cell_value=u_class[1::2], figsize=(15, 12))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To better understand what goes on under the hood of the mpsa class, we manually create the lhs and rhs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of ncasym eliminated:  8\n",
      "max igrad:  2.0\n"
     ]
    }
   ],
   "source": [
    "constit = pp.FourthOrderTensor(mu, lam)\n",
    "stress, bound_stress, *_ = pp.numerics.fv.mpsa.mpsa(g, constit, bound)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "MPSA returns the two sparse matrices \"stress\" and \"bound_stress\". They give define the discretization of the cell-face traction:\n",
    "\\begin{equation}\n",
    "T = \\text{stress} \\cdot u + \\text{bound_stress} \\cdot u_b\n",
    "\\end{equation}\n",
    "Here $u$ is a vector of cell center displacement and has length g.dim $*$ g.num_cells. The vector $u_b$ is the boundary condition values. It is the displacement for Dirichlet boundaries and traction for Neumann boundaries and has length g.dim $*$ g.num_faces.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now ready to set up the linear system of equations and solve it. We define the body forces $f = 0$. Each row in the discretized system is now\n",
    "\\begin{equation}\n",
    "-\\int_{\\Omega_k} f dv = \\int_{\\partial\\Omega_k} T(n)dA = [div \\cdot \\text{stress} \\cdot u + div\\cdot\\text{bound_stress}\\cdot u_b]_k,\n",
    "\\end{equation}\n",
    "We move the known boundary variable $u_b$ over to the right hand side and solve the system:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = np.zeros(g.dim * g.num_cells)\n",
    "\n",
    "div = pp.fvutils.vector_divergence(g)\n",
    "A = div * stress\n",
    "b = -f - div * bound_stress * u_b\n",
    "\n",
    "u = np.linalg.solve(A.A, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives the same displacement as for the Mpsa class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.allclose(u, u_class)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also retrieve the traction on the faces. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/rbe051/anaconda3/envs/porepy/lib/python3.6/site-packages/mpl_toolkits/mplot3d/axes3d.py:738: UserWarning: Attempting to set identical bottom==top results\n",
      "in singular transformations; automatically expanding.\n",
      "bottom=0.0, top=0.0\n",
      "  'bottom=%s, top=%s') % (bottom, top))\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAAKhCAYAAABJgN/wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xl4VOXd//HPLCEhk7AkYSeyCLKI\ngiFsgijiUlxQlAoKKiJ2EVuttnWtoHWtjwvW5Xlsq1b5CVJbl4JSWWyVJUDYZBPZIeyrJJlMlpnz\n+4OHeSBkssDJOXNO3q/rymWSGSbf+XrOmfmc+z73eAzDEAAAAADAPF67CwAAAAAAtyFoAQAAAIDJ\nCFoAAAAAYDKCFgAAAACYjKAFAAAAACYjaAEAAACAyQhaAAAAAGAyghYAAAAAmIygBQAAAAAm89fw\n/katVAEAAAAAzuCpzp0Y0QIAAAAAkxG0AAAAAMBkBC0AAAAAMBlBCwAAAABMRtACAAAAAJMRtAAA\nAADAZAQtAAAAADAZQQsAAAAATEbQAgAAAACTEbQAAAAAwGQELQAAAAAwGUELAAAAAExG0AIAAAAA\nkxG0AAAAAMBkBC0AAAAAMBlBCwAAAABMRtACAAAAAJMRtAAAAADAZAQtAAAAADAZQQsAAAAATEbQ\nAgAAAACTEbQAAAAAwGQELQAAAAAwGUELAAAAAExG0AIAAAAAkxG0AAAAAMBkBC0AAAAAMBlBCwAA\nAABMRtACAAAAAJMRtAAAAADAZAQtAAAAADAZQQsAAAAATEbQAgAAAACTEbQAAAAAwGQELQAAAAAw\nGUELAAAAAExG0AKAOiolJUWbN2+u8LZ3331XAwYMsLgiAADcg6AFALXogw8+UHZ2tlJSUtSiRQsN\nGTJE8+bNO+3H83g82rhx40m/y8/P1/3336+2bdsqEAjorLPO0vDhw7Vo0aJKH6ugoEDt27c/rTpK\nSko0ceJEdezYUYFAQG3bttXYsWO1devW03q82jJx4kSNHj3a7jJq1datW+XxeFRWVmZ3KQCAExC0\nAKCWvPTSS7rvvvv0yCOPaO/evdq+fbvuvvtuffrppzV+rFhvoouLi3XppZdq1apVmj59uo4ePap1\n69Zp5MiR+uKLL2r0WDUxfPhwffbZZ/rggw/0ww8/aOXKlerZs6fmzJlzxo8N8xHCAMAGhmHU5AsA\nUA1HjhwxAoGAMW3atJj3WbRokdG3b1+jYcOGRvPmzY3x48cbxcXF0dslGa+99prRoUMHo23btsZF\nF11kSDKSk5ONQCBgTJ061fjTn/5kNG/e3CgoKKi0nvKPdfx3GzZsMAzDMA4cOGBce+21RmpqqtGr\nVy/jscceM/r371/hY82aNctISkoytm/fHvPv7dy507j22muNxo0bG2effbbx1ltvRW+bMGGCMXz4\ncGPUqFFGSkqK0a1bN2P9+vXGM888YzRp0sRo3bq18a9//St6/4svvth46KGHjF69ehmpqanG0KFD\njYMHDxqGYRhfffWV0apVq5P+dps2bYxZs2YZX3zxhZGQkGD4/X4jEAgY559/vmEYx/7fjB071mje\nvLnRsmVL49FHHzXKysoqfB7BYNC47bbbjEaNGhmdO3c2nn/++ZP+3s6dO40bbrjByMjIMNq2bWtM\nmjQpelsoFDLuvfdeo0WLFkaLFi2Me++91wiFQifV/fzzzxtNmjQxmjdvbnz88cfGjBkzjI4dOxqN\nGzc2nn766ehjhcNh49lnnzXat29vpKWlGT/+8Y+jPcjMzDQkGYFAwAgEAsaCBQuMd955x7jwwguN\n++67z0hLSzMefvhho3Hjxsa3334bfcy9e/ca9evXN/bt2xfz/yMAoELVyk6MaAFALVi4cKFCoZCG\nDRsW8z4+n08vv/yyDhw4oIULF2rOnDl64403TrrPJ598okWLFmnt2rX6+uuvJUkrV65UQUGBRowY\nodmzZ+vKK69UIBCosqYTH6u88ePHKykpSbt379bbb7+tt99+O+bjzJ49W71791ZmZmbM+4wcOVKt\nW7fWrl279NFHH+mRRx7R3Llzo7f/85//1K233qrDhw/rggsu0JVXXqlIJKKdO3fq8ccf109/+tOT\nHu+9997T22+/rd27d8vv9+uXv/xllc/3Rz/6kR555BGNGDFCBQUFWrlypSRpzJgx8vv92rhxo5Yv\nX64vv/xSf/7znyt8jCeeeEJbt27V5s2bNWvWLE2ePDl6WyQS0bXXXqvu3btr586dmjNnjl555RX9\n61//kiQ9/fTTysnJ0YoVK7Ry5UotXrxYTz31VPTf79mzR6FQSDt37tSTTz6pu+66S5MnT9bSpUv1\nzTff6Pe//722bNkiSfrjH/+oTz75RP/5z3+0a9cuNW7cWOPHj5ek6HZx5MgRFRQUqF+/fpKkRYsW\nqX379tq7d69+97vfaeTIkSfVP2XKFA0ePFhNmjSpspcAgNNQ3URmMKIFANU2efJko1mzZjX6Ny+/\n/LJx/fXXR3+WZMyZM+ek++iEUSjDMIzBgwcbDz74YPTn5cuXGw0bNjRSU1ONc845p1qPVVZWZvj9\nfmPdunXR2x5++OGYI1rjxo0zRowYEfN5bN++3fB6vcbRo0ejv3vooYeM22+/3TCMYyNal112WfS2\nzz77zAgEAtFRpaNHjxqSjMOHDxuGcWxE68TnuGbNGiMhIcEoKyurdETr+N8aNWpU9LY9e/YY9erV\nM4LBYPR3H3zwgXHJJZdU+FzatWtnzJw5M/rzn/70p+jfy8nJMTIzM0+6/zPPPGOMGTPGMAzDaN++\nvTFjxozobTNnzjTatGljGMaxEa2kpKRTnnNOTk70/llZWcbHH39sGIZhdO7c2Zg9e3b0tl27dhl+\nv98oLS01tmzZYkgySktLo7e/8847p9R2vN5IJGIYhmH07NnT+PDDDyt83gCASlUrO/ltzHgA4Frp\n6ek6cOCAysrK5PdXfKj9/vvvdf/99ys3N1fBYFBlZWXq2bPnSfepbNTo+N/ZvXt39OcePXroyJEj\nmj17tsaNG1etx9q/f7/KyspOur1NmzaV/s3vv/8+5u27du1SWlqaUlNTT3q83Nzc6M/NmjWLfl+/\nfn1lZGTI5/NFf5aOLdbRqFGjU2pv06aNSktLdeDAgZg1xLJt2zaVlpaqRYsW0d9FIpGYvdm1a9dJ\nt534/bZt27Rr165ojZIUDod10UUXRf/tiX1s06aNdu3aFf05PT39lOdcvi8FBQXRvzVs2DB5vf83\nEcXn82nv3r0xn2v559SnTx8lJyfr3//+t1q0aKGNGzdq6NChMf89AODMMHUQAGpBv379lJiYqE8+\n+STmfX7+85+rc+fO2rBhg44ePapnnnlGhmGcdB+Px1Pp3xk8eLC+/PJLFRYWVllTrMdq0qSJ/H6/\nduzYEf3d9u3bYz7OZZddpsWLFysvL6/C21u2bKlDhw4pPz//pMdr1apVlTXGUr62hIQEZWRkKBAI\nKBgMRm8Lh8Pav39/9OfyzzkzM1OJiYk6cOCAjhw5oiNHjujo0aNas2ZNhX+3RYsWJz3PE+vIzMxU\nu3btoo9z5MgR5efn6/PPP5d0rA/btm07qe6WLVue1vPPzMzUF198cdLfCoVCatWqVcz/rxX9/vbb\nb9fkyZP1/vvva/jw4UpKSjqtegAAVSNoAUAtaNiwoZ588kmNHz9en3zyiYLBoEpLS/XFF1/ot7/9\nraRjy7I3aNBAKSkp+u677/Tmm29W+bjNmjU76bOvbrvtNrVo0ULDhg3T6tWrFQ6HFQqFTho9qorP\n59MNN9ygiRMnKhgMau3atfrrX/8a8/6XXXaZLr/8cg0bNkxLly5VWVmZ8vPz9d///d96++23lZmZ\nqQsvvFAPP/ywQqGQvv32W/3lL385o2XWJ0+erLVr1yoYDOrxxx/X8OHD5fP5dM455ygUCmnGjBkq\nLS3VU089peLi4ui/a9asmbZu3apIJCLpWHC64oor9MADD+jo0aOKRCLatGmT/vOf/1T4d2+66SY9\n++yzOnz4sHbu3KnXXnstelvv3r2Vmpqq559/XkVFRQqHw1q9erWWLFkiSbr55pv11FNPaf/+/Tpw\n4ICefPLJ0+7Bz372Mz366KPR4LZ///7o6pVNmjSR1+uN+ZloJxo9erQ+/vhjTZ48Wbfddttp1QIA\nqB6CFgDUkgceeEAvvfSSnnrqKTVp0kSZmZl67bXXdP3110uS/uu//ksffPCBUlNTddddd2nEiBFV\nPubEiRN1++23q1GjRpo2bZqSkpL01VdfqWvXrrr66qvVoEEDderUSUuWLNG0adOqXetrr72mgoIC\nNW/eXGPGjNEdd9xR6f0/+ugjXXXVVRoxYoQaNmyobt26KTc3V5dddpmkYwstbN26VS1bttSwYcP0\nxBNPRG87HbfeeqvGjBmj5s2bKxQK6dVXX5V0LNC+8cYbGjdunFq1aqVAIKDWrVtH/92Pf/xjScem\n6WVlZUk6trBGSUmJunbtqsaNG2v48OEnTb880eOPP67WrVurXbt2uuyyyzR8+HAlJiZKOhZQp0+f\nrhUrVqhdu3bKyMjQuHHj9MMPP0iSHnvsMWVnZ+v888/Xeeedp6ysLD322GOn9fzvvfdeDR06VFdc\ncYVSU1PVt2/f6OekJScn69FHH1X//v3VqFEj5eTkxHyczMxMZWVlyePxRKc4AgBqh6f8NJUq1OjO\nAACcqUsuuUSjR48+5ZozO7z55puaOnVqzBEwJxg7dqxatmx50gqIAIAaqXxe//9iRAsAgBh2796t\n+fPnKxKJaP369XrxxRcrXbI/3m3dulX/+Mc/dOedd9pdCgC4HkELAIAYSkpK9NOf/lSpqam69NJL\ndd111+nuu++2u6zT8rvf/U7dunXTb37zG7Vr187ucgDA9Zg6CAAAAADVx9RBAAAAALADQQsAAAAA\nTEbQAgAAAACTEbQAAAAAwGQELQAAAAAwGUELAAAAAExG0AIAAAAAkxG0AAAAAMBkBC0AAAAAMBlB\nCwAAAABMRtACAAAAAJMRtAAAAADAZAQtAAAAADAZQQsAAAAATEbQAgAAAACTEbQAAAAAwGQELQAA\nAAAwGUELAAAAAExG0AIAAAAAkxG0AAAAAMBkBC0AAAAAMBlBCwAAAABMRtACAAAAAJMRtAAAAADA\nZAQtAAAAADAZQQsAAAAATEbQAgAAAACTEbQAAAAAwGQELQAAAAAwGUELAAAAAExG0AIAAAAAkxG0\nAAAAAMBkBC0AAAAAMBlBCwAAAABMRtACAAAAAJMRtAAAAADAZAQtAAAAADAZQQsAAAAATEbQAgAA\nAACTEbQAAAAAwGQELQAAAAAwGUELAAAAAExG0AIAAAAAkxG0AAAAAMBkBC0AAAAAMBlBCwAAAABM\nRtACAAAAAJMRtAAAAADAZAQtAAAAADAZQQsAAAAATEbQAgAAAACTEbQAAAAAwGQELQAAAAAwGUEL\nAAAAAExG0AIAAAAAkxG0AAAAAMBkBC0AAAAAMBlBCwAAAABMRtACAAAAAJMRtAAAAADAZAQtAAAA\nADAZQQsAAAAATEbQAgAAAACTEbQAAAAAwGQELQAAAAAwGUELAAAAAExG0AIAAAAAkxG0AAAAAMBk\nBC0AgKssWbJE559/vkKhkAoLC3Xuuedq9erVdpcFAKhjPIZh1OT+NbozAAB2eOyxxxQKhVRUVKTW\nrVvr4YcftrskAIB7eKp1J4IWAMBtSkpK1KtXLyUlJWnBggXy+Xx2lwQAcI9qBS2mDgIAXOfgwYMq\nKChQfn6+QqGQ3eUAAOogRrRQZ4XDYUUiEfn9fnk81ToxAcAhhg4dqpEjR2rLli3avXu3XnvtNbtL\nAgC4R7XeOPpruwogXkUiEW3ZskVt2rRRvXr1CFuAS7z33ntKSEjQLbfconA4rAsvvFBz587VpZde\nandpAIA6hBEt1FmlpaWaN2+e+vbtK4/Ho3r16snrZTYtAAAAKsU1WkB1HA9XxcXFKisrUw1PPgAA\nAACnIGgBkjwejzwej7Zs2aLS0lLCFgAAAM4IQQv4Xx6PR3l5eQqHwyopKVEkErG7JAAAADgUQQso\nx+v1KhKJaNOmTSorK7O7HAAAADgQQQuogNfrVV5enkpLS5lKCAAAgBojaKHOCoVCCgaDKi4ujnkf\nj8ejsrIybdq0ibAFAACAaiNooc5KSkpSvXr1lJubq0OHDlV4n+OLZOzYsUPFxcVctwUAAIBq4QOL\nUaf5/X5lZWXp22+/VVpaWoX3OfGDjIuLi5WQkCCfz8cHHAMAACAmRrRQ5yUlJSk7O1uGYaiwsDDm\nVMLaWgJ+5syZ6tSpkzp06KDnnnvOlMes68aOHaumTZuqW7dudpfiGjt27NCgQYPUtWtXnXvuuZo0\naZLdJblCKBRS79691b17d5177rmaMGGC3SW5Rjgc1gUXXKBrrrnG7lJcpW3btjrvvPPUo0cPZWdn\n210OENcIWoCOLX7RsWNHJSYmKjc3VwcPHqzwfmYvAR8OhzV+/Hh98cUXWrt2raZMmaK1a9ee0WNC\nGjNmjGbOnGl3Ga7i9/v14osvau3atcrJydHrr7/OtmqCxMREzZ07VytXrtSKFSs0c+ZM5eTk2F2W\nK0yaNEldunSxuwxX+uqrr7RixQrl5ubaXQoQ1whawAn8fr969uypzZs3q7i4OOaolVlLwC9evFgd\nOnRQ+/btVa9ePY0cOVKffvrpaT8ejhk4cGDMqaA4PS1atFBWVpYkKTU1VV26dNHOnTttrsr5PB6P\nUlJSJCm6yinTks9cXl6eZsyYoXHjxtldCoA6jKAFlHN8KqEk5ebmxpxKaMYS8Dt37tSaNWuiP7du\n3Zo3rybatm2b3SW40sUXX6zly5erT58+dpfiCuFwWD169FCDBg10+eWX01cT3HffffrDH/4gr9er\nxYsX212Oq3g8Hl1xxRVq2LCh3nrrLbvLAeIaQQuogMfjUWJiotq1a6fc3NxKR63OdAn4UCh0JqWi\nEnzgtPkKCgqUm5urV155RQ0aNLC7HFfw+XxasWKFunXrpsWLF2v16tV2l+Ro06dPV9OmTdWzZ09J\nUklJic0Vucu8efO0bNkytW3bVq+//rq+/vpru0sC4hZBC6hERkaGevbsqeLiYm3YsKHCIHUmS8C3\natVKpaWl0Z/z8vLUqlUrU2oHzFZaWqobb7xRaWlpuuGGG+wux3X8fr8GDRrE9YVnaP78+frss8/U\ntm1bjRw5UkePHtXo0aPtLss1jr9GJSQkaNiwYYwYApUgaAFVSEpKUiAQkHRsKmFFQerEayqOX7dV\nndGtXr16KRQKacuWLSopKdHUqVM1dOhQ84oHTGIYhu6880516dJFzZo1s7sc19i/f7+OHDkiSYpE\nIpo1a5Y6d+5sc1XO9uyzzyovL09bt27V1KlT1aBBA02ePNnuslyhsLBQ+fn5ko5Nef3yyy9Z3RWo\nBJ+jBVRTx44ddfDgQS1dulQHDhxQRkbGKfc5viphnz59dPjw4Wo/dvv27aPf86JVMx6Pp9JQy8IC\ntYO+1p758+eb9vEROKai7bWqYweqZ8iQIVXep3Hjxjp06JAF1QDxhREtoAbS09MVCAS0ZcsWbdiw\nIeb9Dh8+rGAwqA0bNigSicgwjJhfPXr0qPR2q74+/fTTKmt10teRI0fUs2dP2+so//X9999rzZo1\nttdR0y9Jcd3Xir4+/fRT22uo7peT+lr+a926dfruu+9sr8NtfS3/FYlE4mabrmlfa3LiEXATRrSA\nGvJ4PMrOztamTZtUWFioUCikpKSkU+7n9Xq1Y8cO3X333Tpw4EDMxzt06FB02Ww7rVixQueff768\nXuedf6mof3fddZfWrl0bF7090Z49exSJRNSyZUu7S6mxrKysuO1rRZYvX64ePXo4YvStfF+d0N/j\ndu/eLenYRwDEm6q2Vyf12TAMrVixQhdccIHdpVTrOLB9+/ZKX/uAuoCgBZwGj8ejDh06aNeuXVq6\ndKk6deoU877Tp0+vdPW7Cy64QPPmzauNMmskIyNDs2fPVnJyst2l1EggEKiwfzk5OercuXNc9PZE\nL7zwggoKCvTEE0/YXUqNHO9zvPa1Iqmpqfr666/l8/nsLqVKJ/Y11jYdr5566il5vV498sgjdpdy\nisq2V6f1uaysTGlpaXFRc3WOAwMGDLCwIiA+Oe/UNRBH/H6/srOztWXLFoVCoRqtOAgAAAD3ImgB\nZygxMVHZ2dnyeDwxVyWsTMOGDWupMqSnp9tdgivR19pBX2sHfa0d9BWoGkELMMHxDzju0KGDgsFg\njf7tc889V0tV4dNPP7W7BFeir7WDvtYO+lo76CtQNYIWYKK0tLToZ259//33MgyWDgYAAKiLWAwD\nMNnxFc58Pp8KCwu1YsWKSgNXUVGRVq1aZVV5MUUiEa1Zs6bCFRTjXUX9i5e+lrdnz564ra0qq1at\nclztq1atcsRiGOX76qQe79u3Tx6PJy5rrmp7jceaYzm+qFI81Fyd4wAfvA0QtIBac/bZZ2v37t3q\n2LFjpasO5ufnn/SBxXbxer1q27at41YdlFRh/+Klr+WlpaWpsLAwLmurSvv27eO2r7G0b9/eEUGr\nfF+d1ONGjRrJ6/XGZc1Vba/xWHMsx19H4qHmWH0dMWJE9IOJN27cqOzs7OhtJ34vHVvpdubMmbVb\nKGAzghZQi7xerwKBQKVB6/h94kEgEHBk0Kqof/HU1xPVq1dPJSUlcVlbVQKBQNz2NZZAIOCIoFW+\nr07qcb169eJ2u6iqrnisOZbjryPxUHOsvk6fPj36/YABA5SbmytJ0cWiKtO2bVulpqbK5/PJ7/dX\neX/ACQhaAAAAsN1XX32ljIwMu8sATMNiGAAAAABgMoIWAAAAbOXxeHTFFVeoZ8+eeuutt+wuBzAF\nUwcBAABgq3nz5qlVq1bat2+fLr/8cnXu3FkDBw60uyzgjDCiBQAAAFu1atVKktS0aVMNGzZMixcv\ntrki4MwRtAAAAGCbwsJC5efnR7//8ssv1a1bN5urAs4cUwcBAABgm71792rYsGGSji1jf8stt+hH\nP/qRzVUBZ46gBQAAANu0b99eK1eutLsMwHRMHQQAAAAAkxG0AAAAAMBkBC0AAAAAMBlBCwAAAABM\nRtACAAAAAJMRtAAAAADAZAQtAAAAADAZQQsAAAAATEbQAgAAAACTEbQAAAAAwGQELQAAAAAwGUEL\nAAAAAExG0AIAAAAAkxG0AABAnfH222/r8OHDdpcBoA7w210AAAB1XW5urho3bqxOnTrZXYqrTZky\nRb/4xS+0Zs0au0sBUAcwogXAsfbu3Wt3Ca72t7/9Tb169dLLL79sdymu9+qrr6pPnz76zW9+Y3cp\nrjZx4kT5/X796U9/srsUAHUAQQuAo2zevFm//OUv1a9fP3Xt2tXuclytrKxM3333nR577DG7S3G9\nAwcOSJLeeOMNmytxtyVLlmj27Nn6+c9/bncpAOoAghYARzly5IjefvttrV69Wk8//bTd5bjaueee\nq0AgoHfeecfuUlzvgQce0IIFC1RYWGh3Ka7WoEED9erVS88//7zdpQCoAwhaABwlKytLv/71r5WZ\nmalx48bZXY6rnX/++dq1a5duuukmu0txvcGDBzNCCwAuQ9AC4DgTJkzQ0qVL5feznk9t83p5mQAA\n4HTwCgrAcTwej+rXr293GQAAADERtAAAAADAZAQtAAAAADAZQQsAAAAATEbQAgAAAACTEbQAAAAA\nwGQELQAAAAAwGUELAAAAAExG0AIAAAAAkxG0AAAAAMBkBC0AAAAAMBlBCwAAAABMRtACAAAAAJMR\ntAAAAADAZAQtAAAAADAZQQsAAAAATEbQAgAAAACTEbQAAAAAwGQELQAAAAAwGUELAAAAAExG0AIA\nAAAAkxG0AAAAAMBkBC0AAAAAMBlBCwAAAABMRtACAAAAAJMRtFBnlZWVqbi4WKWlpXaXYrvCwkIZ\nhqGCggK7S3E1wzBUWFioYDCoSCRidzmuVlBQwDZtgUgkomAwyDZtAbZpwHn8dhcA2MXj8UiSFi9e\nrPT0dLVp08bmiuyTnZ2tUCik7t27a/fu3XaX41pTpkzRCy+8IEnq2LGjfvKTn9hckXt1795dhmHo\nvPPO0/bt2+0ux7X++te/atKkSZKObdO33367zRW513nnnadIJKIePXpo48aNdpcDoBoY0UKd5fP5\nlJiYqH79+qlRo0b69ttvFQwG9cMPP9hdmuVGjRolSbrxxhttrsTdrrrqKiUkJMjv9+v666+3uxxX\nu+WWWyRJN998s82VuNu1114rv98vv9+va665xu5yXO2mm26SJI0ePdrmSgBUF0ELdZ7X61Xz5s3V\np08fJSYmavPmzVq8eLFKS0tlGIbd5Vnil7/8pRo0aKCHHnrI7lJcrVGjRrr++ut12WWXqWnTpnaX\n42q/+tWvlJqaqgceeMDuUlwtIyNDQ4YM0dVXX6309HS7y3G13/72t0pNTdV9991ndykAqompg8AJ\nfD6fLrjgAgWDQeXk5GjhwoVq3bq16wNXgwYNmDJokXfffdfuEuqEtLQ07dmzx+4y6oSpU6faXUKd\n0KxZM7ZpwGEY0QIqkJycrKSkJGVnZ6usrEyFhYXasGGDiouL7S4NAAAADkDQAipRr149tW/fXikp\nKUpOTtayZctUVFSk/Px8u0sDAABAHGPqIFBNrVq1UsuWLfXNN99o/fr1ko4tEW8YRnQFw+M8Ho+S\nkpLsKLNOCgQCdpdQJ9Bn69Bra9Dn2nX8tbH8ayRQVxC0gBrweDzy+/3Kzs5Wfn6+Fi9erJycHJ11\n1lkn3c8wDIVCIS1YsED9+/dXWVlZzMfMyclR3759a7v0Kl1xxRWaPn266tWrZ3cpNRIIBFRYWHjK\n7+Olr+VNnTpVBw8e1Pjx4+0upUaO9zle+1qRQYMGac6cOfJ643/yxol9jbVNx6u33npL9evX1623\n3mp3KaeobHt1Wp/Lysp05ZVXas6cOXaXUq3jwIABA7Rs2TJJBC3UXfH/6gPEqdTUVNWvX19ZWVkK\nBoMqKCjQpk2bHLtwxsqVK7kGzQJ79uzRjh077C6jTsjNzXXs/ugkO3bsYJEGC5SWlmrFihV2lwGg\nBghawBlKTExUx44dFQgElJA/ZYNOAAAgAElEQVSQED1DGgwGba6sZjwejyKRiN1luJ7P56PPcJVI\nJCKfz2d3Ga4XiUQYGQIchqmDgEk8Ho/OOuss5eXlSZJWr16tYDCob7/9ttKz6kVFRVqzZo1VZcYU\niUS0bt06paam2l1KjVXUv3jpa3n79u3TgQMH4rK2qqxZsyZu+xrLmjVrHBECyvfVST3ev3+/IpFI\nXNZc1fYajzXHEgwGZRhGXNRcnePAOeecY1E1QPwiaAG1pHfv3po3b57atm2rcDgc835Hjhw55Rov\nO/j9frVq1UqNGjWyu5Qaq6h/8dLX8tLS0uK2tqqcddZZjqv9rLPOckTQKt9XJ/U4JSVFaWlpcVlz\nVdtrPNYcS0FBgXw+X1zUHKuvN998sw4ePChJ2rx5s7Kzs6O3nfi9dOzDrmfOnFm7hQI2I2gBtcjr\n9apBgwaVLobh8/niYhTJ6/UqOTk5LmqpqYpqjpe+lpecnCy/3x+XtVUlNTU1bvsay/Ga4135vjqp\nxz6fT/Xr14/LmqvaXuOx5lgikYi8Xm9c1Byrr9OnT49+P2DAAOXm5ko6NuPj+PexhMNhZWdnq1Wr\nVic9DuBkXKMFQNKxoMXCAbXP6/VyjRZcxTAMR6zs6HRu7/OkSZPUpUsXu8sATOXePRZAjbAYhjUI\nWnCbcDjs6gAQL9y8GEZeXp5mzJihcePG2V0KYCqOjAAkHQtajGjVPoIW3IZVB61hGIZrg9Z9992n\nP/zhDwR2uA5bNABJjGhZxev1Vro4CuA0x68dQu1y69TB6dOnq2nTpurZs6fdpQCmc98eC+C0uPVM\nabxhRAtuQ9CyhlunDs6fP1+fffaZ2rZtq5EjR2ru3LkaPXq03WUBpuDICEASi2FYhaAFtyFoWcOt\nUwefffZZ5eXlaevWrZo6daouvfRSTZ482e6yAFNwZAQgiamDViFowW3cOtISb+gz4Dx8jhYASYxo\nWYWgBbdhRMs6bg9al1xyiS655BK7ywBMw5ERgCRGtKxCoIXbsLy7Ndy6GAbgZuyxACQRAKzCiBbc\nhuXdrcHUQcB5CFoAJDGiZRWCFtyGkRZr0GfAedhjAUjiA4utQtCC23CNljUY0QKchyMjAElMHbQK\nH1gMtyFoWcOty7sDbsaREYAkpg5axefz0We4CkHLGkwdBJyHPRaAJEa0rMLUQbgNQcsaTB0EnIcj\nIwBJjGhZhT7DbVje3RqMaAHOwx4LQBIjWlZhRAtuw/Lu1mBEC3AeghaAKAJA7SNowW0IANZgMQzA\neQhaACSxvLtVCFpwG67RsgZTBwHnYY8FIIkAYBWmaMJtCFrWYOQQcB6OjAAkMaJlFT5HC25D0LIG\nUwcB5+HICEASIy1W4XO04DZMabMGI1qA83BkBCCJES2rMEUTbsPy7tZgRAtwHo6MACQxomUVghbc\nhuXdrcHIIeA87LEAJPFBulYhaCGWr7/+Wnl5eXaXUWNOvEbryJEjkqRQKGRzJdXH1EHAeZx1ZARQ\naxjRsgZBC+VFIhG9/vrrGjJkiN555x27y6kxpwWt5557TpmZmZKkJUuW2FxNzTipzwAIWgD+FyNa\n1qDPKO/o0aP697//La/Xq7lz59pdTo05LWjdfPPNSkpKkiT17t3b5mqqjxEtwHmcc2QEUKucuBjG\nqFGj7C6hxpy4vPvKlSvtLsHVGjVqpClTpuj555/X6tWr7S6nxpwWANq0aaPbb79dkpSYmGhzNdXH\nYhiA8xC0AEhy3tTBTZs26Z///KfdZdSYk6YOHjhwQA8++KAGDhxodymu5/f7dffddysnJ8fuUmrM\naSNakvT000/bXUKNsRgG4DzssQAkOWtKm2EYGjFihONGhqRjn6PllED76quv6vXXX3dMvW5w9tln\n211CjTlxeXcnjWQd57SRQwCS3+4CAMQHJ41oRSIRDRw4UImJiVqxYoXd5dSIk0a0JkyYIMMwNHny\nZO3bt8/uchCnWN7dGoxoAc7DHgtAkrNGtHw+n1566SXNnz/f7lJqzElBy+fz6fe//702bdpkdymI\nYwQAazCiBTgPR0YAkpy5GIYTOSloHcebaFTGiddoORGLYQDOw5ERgCSCllWcGLSAyhC0rMHIIeA8\n7LEAJDlr6qCTEbTgNgQtazB1EHAejowAJDlrMQwnc+LnaAGVIWhZgxEtwHnYYwFIYkTLKj6fjz7D\nVQha1uC4ATgPR0YAkrhGyypMHYTbOPFztJyIxTAA5+HICEASUwetwsgh3IbP0bIGUwcB52GPBSCJ\nES2rMKIFt2HqoDVYDANwHo6MACQRtKzCyCHchgBgDaYOAs5D0AIgiQBgFVYdhNswomUNpg4CzsMe\nC0AS1w5ZhVUH4TYEAGswogU4D0dGAJIY0bIK12jBbRjRsgaBFnAe9lgAkhjRsgpBC25D0LIG18IB\nzsOREYAkFsOwCkELbsPnaFmDqYOA83BkBCCJqYNWIWjBbfgcLWswdRBwHvZYAJKYOmgVghbchqmD\n1mDqIOA8HBkBSGJEyyos7w63IWhZgxEtwHnYYwFEMdJS+44HWkIt3IKgZQ1GtADn4cgIQBKLYVjF\n4/HQa7gKQcsaLIYBOA9HRgCSmDpoJa7TgpsQtKzB1EHAedhjAUhiMQwrEbTgFsdPzjDSUvuYOgg4\nD0ELgCRGtKxE0IJb8Bla1mHqIOA8HB0BSGJEy0qsPAi34DO0rMOIFuA8BC0AklgMw0o+n49QC1fg\n+izrMKIFOA9HRwCSxJslCzF1EG5B0LIOi2EAzsMeC0B/+ctf9NFHH+nRRx/VwoUL7S7H1T7//HOV\nlpbqww8/ZASxFi1dulSGYWjJkiV2l+JqH374ocLhsKZPn253Ka6Wm5urBx98UJ988onefPNNu8sB\nUE0ELQA655xzVFRUpKNHj6pDhw52l+NqzzzzjIqKivTEE08QtGrRyJEjZRiGbr75ZrtLcbUJEyao\npKRETz75pN2luFq7du2Un5+vUCikTp062V0OgGoiaAHQRRddpKZNm2rgwIFq0qSJ3eW42sSJEyVJ\n99xzD9OAatGvf/1rSdL9999vcyXu9qtf/UqSCFq1LD09XYMHD1Z6eroGDRpkdzkAqslvdwGAXQzD\nUDgcVllZmfx+doUlS5YoEAjYXYbrDR48WMOHD9f48ePtLsXVxowZo3fffVfjxo2zuxRX+8lPfqLl\ny5fryiuvtLsU15syZYqOHj3KghiAg/DuEnVWWVmZSkpKlJubq3A4rISEBBUVFWnLli0KBAKKRCJ1\n6kLv9PR0u0uoEzwej/7617/aXYbrJSYmcr2hBQKBgN577z27y6gTEhMTmXEAOAxBC3VWQkKC6tev\nr759+0qSSkpKtGjRItWrV09HjhxRKBRSTk6ODMNQUVGRvvvuOwUCAZWVlSkUCikxMdHmZwAAAIB4\nRdAC/le9evXk8/nUqlUrSdKBAwd04YUXyjAMzZ8/XxkZGQoGgyorK9Pq1atVXFysYDColStXKjk5\nWaWlpfrhhx+in3WSlJRk8zOqO5jyaA36bB16bQ36XLuOT3NkuiPqKoIWUAWPxyOv16uMjAxJUl5e\nnrKzsyVJCxYs0Nlnn63CwkLl5eVpx44dCgaDMgxD//nPf1RQUKDmzZtXurrc3r171axZM0ueS2UW\nLVqk3r17O+4FMTs7W7m5uaf8Pl76Wt6ePXsUCoXUtm1bu0upkeN9jte+VmTRokXq06eP3WVUy4l9\njbVNx6vt27fL7/erZcuWdpdyisq2V6f1WYqfbbo6x4Hnn39e06ZNk0TQQt1F0ALOUEpKilJSUrRp\n0yZ169ZNR48elSRlZWUpJydHrVu3VjgcjvnvDx06FBdvUH77299q2bJlSk1NtbuUGquof/HS1/I+\n/fRT7du3TxMmTLC7lBpr2bJl3Pa1Ir/4xS+0bds2+Xw+u0upUvm+OqXHkvTOO++oQYMG+sUvfmF3\nKaeoant1Up+Lior0wAMPaPPmzXaXErOvt9xyiw4dOiRJ2rx5c/SkpKSTvpekjIwMzZw5s3YLBWxG\n0AJqSUJCgnw+nxo1aqSysrKY9/P5fGrYsKGFlcWuIzU1NS5qqamKao6XvpaXnJwsv98fl7VVpWHD\nhnHb11iO1xzvyvfVST32+/1KTk6Oy5qr2l7jseZY/H6/vF5vXNQcq68zZsyIfj9gwIDoiKHH46l0\n9DAUCmngwIEqLi5WWVmZhg8frieeeML8wgGLEbQASJK8Xq8ikYjdZbiez+ejz3CVSCTiiDDrdG5e\nBTcxMVFz585VSkqKSktLNWDAAA0ZMiS6WBXgVO7cYwHUGEHLGvQZbhMOh10bAOKJm4OWx+NRSkqK\nJKm0tFSlpaVc1wVXcOceC6DGPB4PAcACjGjBbQzDYETLApFIxNXhIxwOq0ePHmratKkuv/zyuFj0\nAzhTBC0Ako4FrcpWR4Q5vF5vpYujAE7DiJY1DMNwdZ99Pp9WrFihvLw8LV68WKtXr7a7JOCMuXeP\nBVAjTGmzhs/nI2jBVcLhMCNaFnDz1METNWrUSIMGDWJFQriC+/dYANXi9XoZ0bIAgRZuw2IY1nDz\n1MH9+/fryJEjko4tYz9r1ix17tzZ5qqAM8eqgwAkEQCswtRBuA1TB63h5qC1e/du3X777QqHw4pE\nIrrpppt0zTXX2F0WcMYIWgAkEbSswtRBuA1Byxpuvkbr/PPP1/Lly+0uAzCdO/dYADXGqoPWYNVB\nuA1TB63h5qAFuBV7LABJXKNlFUYO4TaMaFmjriyGAbgJeywASQQAqzCiBbdhRMsaBC3AedhjAUgi\naFmFxTDgNgQAa7h5MQzArTgyApBE0LIKi2HAbfgcLWtwjRbgPOyxACSxGIZVmDoIt2HqoDUYOQSc\nhz0WgCQWw7AKUwfhNiyGYQ2CFuA87LEAJDGiZRWmaMJtGNGyBtdoAc5D0AIgiQBgFUa04DbhcJgA\nYAFGtADnYY8FIImgZRUWw4DbsBiGNVgMA3Ae9lgAkpg6aBWfz8e1cHAVpg5ag6mDgPMQtABIYjEM\nqzB1EG7DlDZr0GfAedhjAUgiaFmFqYNwG6YOWoMRLcB5CFoAJHGNllX4HC24DVMHrcE1WoDzsMcC\nkETQsgpTB+E2fI6WNQhagPOwxwKQRNCyCiNacBtGtKzBNVqA87DHApDEqoNW8Xg8jGjBVRjRsgZB\nC3Ae9lgAklgMwyoshgG3IWhZg8UwAOfhyAhAElMHrcLnaMFtmDpoDa7RApyHPRaAJKYOWoXFMOA2\njGhZg6mDgPOwxwKQxIiWVZg6CLdx8ojWtGnTHDPCTNACnIc9FoAkRrSsQqBFdbz00ksqLi62u4xq\ncXIAuOuuu3THHXfYXUa1cI3WqR5//HG98sor0Z8fffRRTZo0ycaKgJM588gIwHRODAD79u2zu4Qa\nc+KIVjAYlCS98MILeuaZZxwzAuBEs2fPliRNmDBBS5cutbma6gmHw44b0ToeYlNTUzV8+HCbq6ke\nwzAIWuWMHTtW7733nqRjQXTq1KkaPXq0zVUB/8dvdwEA4oMTVx3Mzs6WJO3fv19NmjSxuZrqcWLQ\n+uabbyRJTzzxhFJSUnizV4uuu+46ScdGmJctW6YLL7zQ5oqq5sQRre+++y7635SUFJurqR4n9rm2\ntW3bVunp6Vq+fLn27t2rCy64QOnp6XaXBUSxxwKQ5Myg9fDDD0uSunTpotzcXJurqR4njhxeeeWV\nkqRWrVqpc+fONlfjboWFhZKkGTNmaNiwYTZXUz1OHNHq3r27JDkmZEkErVjGjRund999V++8847G\njh1rdznASRjRAiDJmQHgrLPOkiT17t1bbdq0sbma6nHyqoPLli3TvHnz7C6jTrjooovsLqHanBi0\nnIigVbFhw4bp8ccfV2lpqT744AO7ywFOQtACIMmZQWvQoEGSjp39d8p0Np/P57g+HxcIBNSwYUO7\ny0CccfKqg05C0KpYvXr1NGjQIDVq1IjtEHGHoAVAkjNXHUxOTpYkx4QsyZmBFqhMOBx21D7oVCyG\nUbFIJKKcnBz97W9/s7sU4BScGgEgyZnXaDmRExfDACpjGAYjCRYwDIMRrXLWrl2rDh06aPDgwerY\nsaPd5QCnYEQLgCRGWqzi5KmDQEXC4TABwAJMHTxV165dtXnzZrvLAGJijwUgiaBlFScvhgFUhMUw\nrEHQApyHPRaAJIKWVRjRgtuwGIY1CFqA87DHAogiANQ+RrTgNkwdtEYkEmExDMBhODICkMRiGFZh\nMQy4DSNa1mAxDMB52GMBSGLqoFWYOgi3YUTLGkwdBJyHPRaAJIKWVY5P/aHXcAsWw7AGUwcB5yFo\nAZBE0LIS0wfhFoZhMKXNIoxoAc7DHgtA0rGRFq7RsgbTB+EWx0dZGGmpfYxoAc5D0AIgiREtK7Hy\nINyChTCsw4gW4DzssQAkseqglZg6CLdgIQzrMEUTcB72WACSGNGyEr2GW7AQhnUIWoDzsMcCkMSb\nfysxogW3YOqgdZg6CDgPeywASQQtK7EYBtyCqYPWIWgBzsMeC0DSsVUHefNvDY/Hw4gWXIERLeuw\n6iDgPAQtAJJY3t1KjGjBLRjRsg7XaAHOwx4LQBJTB61E0IJbELSsw9RBwHnYYwFIImhZic/RglsY\nhsHUQYsQtADnYY8FIImgZSVWHYRbMKJlHa7RApyHoyMASSyGYSWmDsItWAzDOgQtwHkIWgAkHRvR\nYjEMazB1EG7BiJZ1WAwDcB72WACSCFpWYkQLbhEOhxnRsgjXaAHOwx4LQBLXaFmJES24BVMHrUPQ\nApzHb3cBAOz3//7f/9MLL7wgSWrTpo1Gjhxpc0XuRtCCWzB10Bp///vfNWHCBElShw4dNGbMGHsL\nAlAtHB0B6JJLLpHX65XX69XFF19sdzmu9uGHH2rPnj2aNGmSCgoK7C7Htfbv3y/DMLR//367S3Gt\nYDCol19+Wfv27dMHH3xgdzmuNmDAAPl8Pnm9Xg0aNMjucgBUE0ELgFq1aqV+/fqpd+/eatGihd3l\nuNpXX32lgwcP6uOPP7a7FFfr06ePDMNQr1697C7F1f7xj3/o8OHDmjVrlt2luFqzZs3Uv39/ZWdn\nq02bNnaXA6CaCFoAJEkzZ87kzZIFHnvsMUnSTTfdpJSUFJurca+xY8dKku644w6bK3Gv5ORk3XLL\nLZKkiRMn2ltMHTBjxgzNmTPH7jIA1ADXaAEmMgwjuqDE9u3bFQqFtG7dukoXmQiFQlq/fr1VJbpS\nRf2L577ef//9uuqqq+K2vljWr18f13090ZAhQ/Tmm2/qmmuucUS95fvqhJolafTo0UpNTY3b7aKq\nuuKxZieozv/v9u3bW1QNEL8IWsBpKC0tVWFhoUpLS7VhwwYVFhaqoKBACxcuVCgUknTsA4D9fr+a\nNWtW6cIH+/fvV9OmTa0q3ZUq6l889/X++++3u4TT0rRp07ju64maNm2q1atX211GtZXvqxN6LB2r\n89e//rXdZcRU1fbqlD7Hm1h9HTVqlA4fPixJ2rx5s7Kzs6O3nfi9JGVkZGjmzJmSpB07dui2227T\n3r175fF49JOf/ET33ntvLT4DwBoELaAS+fn5CgaDKi4u1qpVq6KBavny5QoEAopEImrYsKFatmyp\nYDCoCy+8UAsWLJAkZWZmaseOHUpLS1NZWVnMv+H3+9W4cWOrnpIrVdQ/+mq+xo0b09daUr6v9Ngc\nVW2v9Pn0xOrr559/Hv1+wIABys3NlXTsxOPx72M93osvvqisrCzl5+erZ8+euvzyy9W1a1fziwcs\nRNBCnWYYhoqKilRYWKjCwkKFQiHl5uaquLhYwWBQmzdvViAQkNfrVZs2bZScnKzFixerd+/ekqQf\nfvihwrN6Ho9HSUlJVj+dOisQCNhdQp1An61Dr61Bn2uXx+M56b+xtGjRIroQU2pqqrp06aKdO3cS\ntOB4BC3UWcen/61bt06BQEDJycny+/0677zzVK9ePS1cuFDdu3eXJO3bt08NGjSo9mMbhqFQKKQF\nCxaof//+lY5o5eTkqG/fvmf8fM7U0KFD9fe//10JCQl2l1IjgUBAhYWFp/w+Xvpa3rRp01RYWOi4\nRRqO9zle+1qRq6++Wv/85z8d8TlPJ/Y11jYdr95//315vV6NGjXK7lJOUdn26rQ+h8NhXXfddZo+\nfbrdpVTrODBgwAAtW7ZMUtVB60Rbt27V8uXL1adPnzOqEYgHBC3UWQkJCUpJSVFWVlb0dzt27FBi\nYqKNVdln4cKFKikpcVzQcppt27YpPz/f7jLqhK+//lqGYdhdhutt3bpVfj9vJ2pbOBzWvHnz7C6j\nVhUUFOjGG2/UK6+8UqOTm0C8iv/TfAAAAHWc208alJaW6sYbb9SoUaN0ww032F0OYAqCFoAot7+Q\nxwv6DDdhe7ZOTabgOYlhGLrzzjvVpUsXx67KClSEoAVAkntfwOMNfYYbsV3jTMyfP1/vv/++5s6d\nqx49eqhHjx4nrWAIOBWTqgFEcWYaQE1x3LCGm/s8YMAAVz8/1F2MaAGQxBlpK/GGAm7D8cMa9Blw\nFoIWAFiIN0oATgcnaADnIWgBiOKF3Br0GW7C9mwdTtQAzkLQAiCJF3Cr0Ge4Edt17SPQAs5D0AIQ\nxQs5gJriuGEdAi3gLAQtAJJ4AbcSb0zhNhw/AOBUBC0AUQSA2scbUgCng+Mz4DwELQCSCAAATg8B\nwDocpwFnIWgBgMV4YwqgpjhuAM5D0AIQxQt57eOMNIDTxfEDcBaCFgBJvIADOD2coLEGfQach6AF\nIIoXcmvQZ7gNJ2qsQZ8BZyFoAZDEC7hV6DPchhMH1qDPgPMQtAAAwBnhBII16DPgLAQtAFGcMbUG\nfQYAwP0IWgAkcabUKvQZbsOJA2vQZ8B5CFoAonghtwZ9httwAsEa9BlwFoIWAEm8gFuFPsNtOHFg\nDfoMOA9BCwAAnBFOIFiDPgPOQtACEMUZU2vQZwA1xXEDcB6CFgBJnCm1Cn2G2xAArMPxA3AWghaA\nKN4wATgdBAAAOBVBCwAsRqCFm7A9W4M+A85D0AIAC3HmH27Edm0N+gw4C0ELQBRnTAEgPnF8BpyH\noAVAEmdKAZweAoB1OE4DzkLQAhDltDdMH330kd0l1JgT3yh99dVXdpdQJ2zfvl1PP/203WWcFqds\n14ZhaNWqVdqzZ48kaenSpTZXVH1OOz4DkPx2FwAgPjjljZIklZaW6sUXX9Rzzz0nSTp48KDS09Nt\nrqr6nPKGaceOHbr66qu1fft2u0txvcGDB0uSvF7nnf90yvYsSYWFhbr11lu1YcMGSdJDDz2kWbNm\n2VxV9TnpOA2AES0ADvQ///M/+vOf/6y0tDRJ0sqVK22uqPqc9EYpPz9fxcXF8vl8dpfiem+88YaG\nDh2qSCRidymnxSnbdUpKipYvX64777xTkjR8+HCbKwLgZgQtAFFOOTN9zz33aOPGjdq8ebMk6dJL\nL7W5Infq2rWr1q9fz4iWBTp16qQpU6bom2++sbsU1/N4PJo0aZIkacSIETZXU31OOT4D+D8ELQCS\nnHNG2g2c9oYpEAjYXUKdkZWVZXcJNea07Vn6v+Ndo0aNbK6kZjhOA85C0AIAC/FGCW7Edl37nBho\ngbqOoAVAEm+UAAAAzETQAhDFGVNr0Ge4CduzNegz4DwELQCSGNGyCn2GG7FdW4M+A85C0AIQxRlT\na9BnADXFcQNwHoIWAEmcKbUKfQZwujh+AM5C0AIAAGeEAAAApyJoAYhiaoo16DPchO3ZGvQZcB6C\nFgBJnJG2Cn2GG7FdW4M+A85C0AIQxRlTa9BnAADcj6AFQBJnSq1Cn+E2nDiwDscPwFkIWgCieMME\n4HQQAGofx2fAeQhaACTxRslKvGECcDo4TgPOQtACEEUAqH28UYLbcNywBn0GnIegBUASAcBKvGGC\n23D8sAZ9BpyFoAUAFuKNEgAAdQNBC0AUIy0AaorjhjXoM+A8BC0AkhhpsRJvmOA2HD+sQZ8BZyFo\nAYgiANQ+3ijBbThuWIM+A85D0AIgiQBgJd4wwW04fliDPgPOQtACEEUAqH28UQJwOjg+A85D0AIg\niQAA4PQQAKzDcRpwFoIWgCjeMFmDPsNtCAC1j+MG4DwELQCSeKNkFfoM4HRx/ACchaAFQGVlZYpE\nIiopKbG7lDqBM9NwE7Zna5SUlCgSiaisrMzuUgBUE0ELgK688kpt3bpVffr0UX5+vt3lAHAYRlpq\nVzAYVHZ2tvLy8jR48GC7ywFQTX67CwDcKhgMqri4WOvXr1ckEol5v+LiYm3YsMHCyk41ePBg5eTk\nqFu3btqzZ4/27Nljaz01VVH/4qGvFdm3b5+OHj0al7VVZcOGDXHb11g2bNggn89ndxlVKt9XJ/X4\nhx9+0L59++Ky5qq213isOZbu3btr6dKlGjx4sO11V+c40K5dO4uqAeIXQQswWTgcliStWrVKHo9H\n6enp0d9VZO/evUpLS7OqvArdcccdevXVV/Xggw/aXsvpqKjmeOhreRs3btS0adN06NAh5eTk6Kqr\nrrK7pBpJS0uLy75WJC8vT4ZhKBgMqk2bNnaXU6XyfXVCjyXpyy+/1MKFC7VmzRpdfPHF6tSpk90l\nnaSq7dUpfZakhx9+WGPGjNG4ceOUmJhoay2x+nrrrbfq0KFDkqStW7cqOzs7etuJ30tSRkaGZs6c\nKUkaO3aspk+frqZNm2r16tW1WDlgLYIWYALDMFRaWqpFixapuLhYktSnTx8tWLBAGRkZlc6p9/v9\nSk9Pt6rUmJw2inWiivoXL3090dGjR7Vu3TpJUkJCQtzVV5X09PS47GtFevXqJcMwdN1112nHjh12\nl1Ol8n11Qo+lY9vx7t27tXv3bjVp0iTu6q5qe423eiszZMgQ7d271+4yJMXu6+effx79fsCAAcrN\nzZV0bGrp8e8rMmbMGN1zzz267bbbzC8WsBHXaAFnwDAMbdu2TQsWLFA4HNZ5552n5ORku8tCnGrX\nrp2ys7NVv359/fjHP7uUnloAACAASURBVLa7HFe7++67JUk/+9nPbK7E3W644QYFAgF1795dHTt2\ntLscONTAgQMdNboIVBcjWsBpCIVC2r59uwoLCxUOh9W7d28tWbKEkIUqvf/++9qyZYsjrhtysp//\n/Od66623dM8999hdiqt5vV794x//UMuWLe0uBQDiDkELqIH8/HwVFRVp+fLlOuuss5SSkqL27dvb\nXRYcJDMzU5mZmXaX4XqBQEDff/+93WXUCf3797e7BACISwQtoAqGYaisrEy5ubnyeDxKSEhQ3759\n5fF4tG3btgr/jcfjUVJSksWV1l2BQMDuEuoE+mwdem0N+ly7ji/7z/L/qKsIWkAl8vLytGPHDpWW\nlqpz585KSUnRggULqnzRMAxDoVBICxYsUP/+/StdDCMnJ0d9+/Y1u/QaW7hwofr162d3GTUWCARU\nWFh4yu/jpa/lbd26VYmJiWrRooXdpdTI8T7Ha18rsmjRIvXp08fuMqrlxL7G2qbj1d69e1VYWBiX\no/uVba9O67MUP8fp6hwHBgwYoGXLlkkiaKHuYjEMoJySkhJt2rRJBQUFKioqUlZWlurXr6+UlBS7\nS6s1oVBI1157rd1l1AmvvfaaPv74Y7vLqBMuv/zySj9aAeb47LPP9Mc//tHuMuqEYcOG6ejRo3aX\nYbqbb75Z/fr10/r169W6dWv95S9/sbskwBSMaAH/q7CwUKFQSLm5uWrdurUCgUCdWUUrFArZ/rks\ndUVRURHTSi1iGIbdJdQJiYn/v717j46yPNc/fs0kmZmEQ0jCQQynhCSEcBDEsDkIbpWq0BZFLILV\nHjzUtu6ibdk1asHKpurettZ2Wd111dXa1oLVWtGqWNBStVgCiiA5E0KAEAI5kEBmJnP8/cFOfooE\nQYf3nXfe72etrgVmFtw8nUzmmud+7sctr9drdhm24Ha7e64QSSSrV682uwTgrCBowdai0aja2tq0\nZ88eBYNBJScna/r06XI4HNq/f7/Z5RnG7/fL5XKZXYYteL1epaamml0GEDNpaWny+Xxml2ELLpcr\nIYMWkKhoHYRtBQIBeb1e7d27V7m5uZo6daqSk5Nt2Uve1dXFLotB/H4/1wAgoaSlpbGjZRCPxyO/\n3292GQBOEztasK2UlBSlpqbqvPPOM7sU0xG0jOP1ellrJBSPx8OOlkEStXUQSFTsaMG2HA6HnE6+\nBSRaB43k8/nY0UJCoXXQOAQtwFp4lwmAHS0DEbSQaFJTUwlaBqF1ELAWghYAdXV1MXXQIEwdRKIh\naBmHHS3AWghaABjvbiCv18uOFhIKwzCMQ9ACrIWgBUCBQICgZRB2tJBoaGczjsfjIWgBFkLQAiC/\n38+bf4NwRguJhh0t47hcLkItYCEELQC0DhqIoIVE072jFY1GzS4l4bF7CFgLQQsArYMGCQaDikaj\nSklJMbsUIGacTqfcbjcBwAAej0eBQMDsMgCcJoIWAFoHDcIgDCQq2geNQesgYC0ELQDq6uriwmID\nMAgDicrj8TDi3QAMwwCshaAFgAuLDcL5LCSqtLQ0gpYBGO8OWAtBCwCtgwYhaCFRpaam0jpoAM7C\nAdZC0AJA66BBaB1EokpNTSUAGIDWQcBaCFoAaB00CMMwkKjY0TIGrYOAtRC0AKirq4vx7gZgRwuJ\nKjU1lTNaBiBoAdZC0ALAhcUG4YwWEhXDMIzBGS3AWghaANjRMojP51NqaqrZZQAxR+ugMTijBVgL\nQQsAZ7QMQtBCoqJ10BgELcBaCFoAaB00CMMwkKgIWsZwuVy0DgIWQtACQOugQRiGgURF0DIGO1qA\ntRC0ANA6aBCGYSBRpaWlcUbLAAQtwFoIWgDk9/u5sNgAnNFComJHyxi0DgLWQtACoEAgwI6WAQha\nSFTxGLRaWlr0/PPPq7a2VpFIxOxyYsLj8SgQCJhdBoDTRNACIL/fT9AygNfrJWghIcXjePeqqir9\n8pe/1HnnnadZs2aZXU5McI8WYC0ELQBMHTSI3+/njBYSUlpaWtwFgBkzZuj8889XNBqV1+vV/v37\nzS7pM3O73ZzRAiwk2ewCAJgvEAgQtAzAjpa9BINBPfnkk/r73/8ul8ul3//+92aXdNZ4PJ6429GS\npO9+97tqampSZ2enLrzwQj3++OOaO3eu2WV9agzDAKyFoAWAHS2DcEbLXjZv3qxly5YpGo3K4/Ho\nuuuuU2FhocaOHavCwkIVFBQkzPddPO5oSdK5556rq6++WpI0c+ZMfe1rX9Obb76p+++/3+TKPh1a\nBwFroXUQAPdoGYSgZS8XXnihtm7dqpycHF1zzTW6+uqr5XA49Je//EVf/epXNXToUE2aNEmLFy/W\nb37zG7W3t5td8qcWj2e0TjR9+nRt2rRJDofDsrtCLpdLoVAoYYZ7AImOHS3A5kKhkKLRqFJSUswu\nJeHFU+tgJBLRO++8o5kzZ/b8N5/Pp+rqalVUVKiyslKVlZUaM2aM7rvvPhMrtbbCwkLt2LFDkUhE\nyckf/ZEbCAS0a9cuVVRUaOPGjQoGgyZV+dlZIWhJUlZWlmV3syTJ4XDI5XKpq6srbl5LAPSOoAXY\nnNUmDnZ1dWn58uXKzMxUSUmJ2eWckXgZhrFnzx4tWrRIZWVluvXWW7V//35VVlaqoaFBOTk5Gjt2\nrMaOHatFixZp2rRpZpdreU6nU07nxxtIXC6XioqKVFRUpOzsbA0cONCE6mIjXlsHE5HH45Hf7ydo\nARZA0AJsrquryzJBa/fu3frKV76iYcOG6a677jK7nDMWLztaCxYsUHV1tRwOhzo6OnTddddpzJgx\nysvLY2cTn4pVdrQSAZMHAesgaAE219XVJZfLZXYZvYpGo1qyZIlyc3P19NNPq6SkRN/85jflcDjM\nLu2MxcsZrXfffVfbtm3TK6+8ouuuu06jR482uyRYXGpqKjtaBiFoAdZB0AJsLt5bB//2t7/pr3/9\nqyTp2WeftfRo5ngJWk6nU1OmTNGUKVPMLgUJgh0t43S3DgKIf0wdBGwu3lsHly1bJknKzc1VZmam\nydV8etFoNG6CFhBrLpdL4XBYoVDI7FISHjtagHWwowXYnN/vj+vWwc2bN8vtdispKcnsUj6T7hZN\nq/87gJNxOBxKS0uTz+dTv379zC4noRG0AOtgRwuwuUAgENc7WmlpaQkRTuJlEAZwtng8HtoHDeDx\neAhagEUQtACb8/v9XFZsANoGkegY8W4Ml8vFOgMWQdACbI6gZQyfzxcXd2gBZ0taWho7WgZgRwuw\nDoIWYHPx3jqYKLxeL+uMhObxeOTz+cwuI+ERtADrIGgBNseOljHY0UKi6x6GgbOL1kHAOghagM0R\ntIzBGS0kOoZhGIMdLcA6CFqAzdE6aAyCFhIdwzCMQdACrIOgBdgcO1rGoHUQiS41NZUdLQPQOghY\nB0ELsLmuri6ClgEYhoFEl5qayhktA3g8HgUCAbPLAHAaCFqAzRG0jMGOFhIdwzCM4Xa72dECLIKg\nBdic3+9np8UAPp+PdUZCYxiGMQhagHUQtACbY0fLGOxoIdExDMMYtA4C1kHQAmyOoGUMghYSHcMw\njMGOFmAdBC3A5rq6umhpMwDDMJDoUlNTCQAGcLvdjHcHLIKgBdic3++Xy+Uyu4yEx44WEl1aWho7\nWgYgaAHWQdACbI4dLWNwYTESHa2DxvB4POwcAhZB0AJsjqBlDIIWEh2tg8ZgRwuwDoIWYHO0DhrD\n6/XSOoiExo6WMQhagHUQtACbY0fLGNxXhkTHjpYxaB0ErCPZ7AKARBUOh9XV1aWamhpFIpFeH9fV\n1aXa2loDK/uojo4ONTc3m1rDZ3Wy2s1e1xO1traqtbU1rmo6U7W1tXG3rp+ktrZWSUlJZpfxiU5c\nVyutcbfW1la1tbXFVe2f9HyNp1pP1+HDh3X06FFTaz+d14GRI0caVA0QvwhawFlSWloqSerfv/8p\ng1ZSUpL69+9vVFkfEwqFlJGRYWoNn9XJajd7XU8UCoU0cODAuKrpTPXv3z/u1vWTdNcc705cVyut\ncbeBAwcqEAjEVe2f9HyNp1pPV2ZmpoLBoKm197auX/nKV9TW1iZJ2rNnjy644IKer33419Lx58u6\ndet6fr9u3TrdfvvtCofDuvnmm1VSUnKWqgeMQ9ACYiwYDEqSioqKVFZWpiFDhigUCvX6+NraWg0a\nNMio8j4mHA5r6NChptbwWZ2sdrPX9USBQEDnnntuXNV0pgYNGhR36/pJBg0aZImgdeK6WmmNu3m9\nXgWDwbiq/ZOer/FU6+k6duyYwuGwqbX3tq6vvvpqz68vvPBCbd26VZLkcDh6fn0y4XBYt912m9av\nX69hw4apuLhY8+fPV1FRUeyLBwzEGS0gRqLRqMrKynqCVnp6uskVfbLrr79eu3fv1le/+lV6/s8y\nhmEg0TEM4+wLBoO6/vrrtW/fPi1evNjscmKmtLRUeXl5ys3Nlcvl0uLFi7V27VqzywI+M4IWEAOd\nnZ3yer3q27evpd5MDx48WOFwWIFAQG632+xyEtb69evV0dGhbdu2mV2KbXzve99TJBLRHXfcYXYp\ntrF9+3YdPXpUr732mtmlJKzk5GRFIhHTd7RiraGhQcOHD+/5/bBhw9TQ0GBiRUBsELSAz6ixsVHb\nt2+Xx+Ox3OHf73//+5KkBx98UA6Hw+RqEteyZcvk9Xr1zW9+0+xSbKM71BJujXPLLbeoq6tL3/3u\nd80uJWE5HA795Cc/kXT8dQVAfCNoAZ9SOByWz+dTU1OTiouLLXEO5ETZ2dn64IMP9LnPfc7sUhLa\n7bffLklasWKFyZXYx49//GNJ0qpVq0yuxD5+9KMfSZK+853vmFtIgps5c6Y++OADy32wdyrZ2dna\nt29fz+/379+v7OxsEysCYoNhGMCn0NnZqR07digpKUnnnXder7tB0WhUXV1dvU4dDIfDikQiPee6\nzDB8+HBT//5YOLH+aDRq+rp+2DXXXKM1a9bo+uuvj5uaPo1gMKhoNGqJf0NxcbF+9atfaebMmZao\n98R1tULNJ1q8eLHWrFmjxYsXx039kUhEgUCg19foeKnzTJn9uh2JROT3+3v9gDESiaipqUmdnZ3q\n06fPJ/55xcXFqqmpUV1dnbKzs7VmzRr98Y9/jHXZgOEIWsAZCgaD2r59u8aPH6+dO3f2+gM8EonI\n6XT2jHk/2df9fr9SUlJOOY0Jn+zE9QsEAopGo3G1rvfdd5/ef/99s8v4TLZs2SKfzxdX63oqo0eP\ntkytPp9PW7Zs6Xk9sUrdJ7r33ntVWVlpdhk9AoGA3nnnHblcrpN+3arrbLZQKKR//vOf8ng8cjr/\nf3NUSUmJOjo6FI1GlZSUpFGjRumcc86RdOrx7snJyXr00Ud1+eWXKxwO68Ybb9S4ceOM+wcBZwlB\nCzhN4XBYFRUVCgaDmjlzppKTT/3tEwqFNGXKlJN+7ciRI6qurtakSZM0YMCAs1GurUydOrXn1+3t\n7dq1a5cmT578kTcA+Oxyc3PV0dGh0aNHm11Kwtm9e7f69u2rwYMHS/rocxqfXiQS0fvvv6/c3NyT\nvtayzp9ee3u7qqqqlJ+fr4yMDEnSG2+88ZHHbNmypad1+sMfJJzMvHnzNG/evLNXMGACghZwGiKR\niEpLS5Wdna2Ojo5Thiyn06l+/fpZ8iJMq2KtjUHAMg7PaWOwzsZwOp0aMWKEhgwZctKvn3h5MZAo\nCFrAJzhw4IC8Xq+mTZum/v37a//+/b0+NhgM6vzzzz/pWNpIJKKamhqFQiEVFhbG1fCMnTt36o9/\n/KPuv/9+s0v51KLRqHbu3KkhQ4b07ArEox/84Ae6+eabVVBQYHYpZ2zr1q2aNGnSJ+7mxotrr71W\nq1evtsTOZjgc1nvvvafi4mKzSzljtbW1evTRR/Wzn/3M7FJ6dfjwYTU2NmrChAmWnrC6YsUKLViw\nQJMnTza7lB7hcFhVVVVyOp0qKCg46fdbKBTSqlWrtHPnTj399NMaOHCgCZUCxrPGT0vABN1v3MPh\nsPr06dPrJ5/RaFTRaPSU57Gi0ah8Pp+Sk5Plcrn03nvvnc3Sz9hLL72k2tpabdmyxexSPrVAINBz\n7q2+vt7scnr14osvatq0aWpvbze7lDPS/Ry20rj01157TaWlpXH1ocapeL1elZaWWi4IHDhwQC+9\n9JKuu+46s0s5Jb/fr02bNvV6XssKamtr9dxzzykUCpldyscEAgG99dZbSk1N/dhzuKSkRO3t7QqH\nwyosLNSgQYN6HZLB7hYSCUELOIljx47J6/VqxIgRGjZsmN55551eH9s9OfD8888/6dePHj2qiooK\njR8/XllZWWer5M/k+eef14wZMyz5abp0fI2rqqo0efLkuH5THQ6H1dzcrHnz5snj8Zhdzhlpbm7W\nkSNHlJeXZ3YpZ8RKVy/U1taqf//+lruINhgM6qabbtJ5550X1yEmEonovffeU0FBgWVbBmfMmKG6\nurq4fa1ubW3Vrl27NGbMmI+s8euvv97z6127dunmm2/WN77xDd1yyy2W+2ABOBPx308BGKyhoUE7\nduyQx+PR8OHDTzm6vXuMeG+amppUWVkZ1yFLkmpqaiz3BrpbKBRSZWWlxo4dG/dvqBsaGjRw4EDL\nhSzp+AAXBrecXRkZGTpy5IjZZZyxlJQUDR069CP3IMUjp9OpoqIiVVZWxuWO0OnIy8tTTU2N2WX0\nKjMzUxMmTFBVVZUOHjx40sfk5eXp5Zdf1uuvv64bb7xRPp/P4CoB47CjBfyf7guIm5ubNXXq1F7b\nAKXjn4yGQiElJSX12m7XfX+Wx+NRWVnZ2So7Jnbu3Cm/32/J1sHulszy8nKzS/lE27ZtU1ZWliXX\n2ev1qrW1VXV1dWaXcka2bNkS9wG8WzQaldfrtWTYyszM1N/+9rePjfCOR6FQSJs2bVJqaqrZpZwx\nv9+vsrKyuH8NiUajqq6u1q5du+R2uz/yte42Qun4v2fQoEHKycn52OMk2ghhfQQtQMdbBT/44AMl\nJSVp4sSJn7iLFQwG5XA4TtouGAqFVFZWpqysLOXk5MR9W0QgEFBLS4u++MUvxnXbz8k0Njaqra1N\nRUVFZpdyWiorKzVhwoS4bfvpTSgU0vvvv2+JN9EnslLroCS9++67mjhxolJSUswu5YxMnDhRbrfb\nMs/tiooKpaen69xzzzW7lDMSDod18803a/z48XEfFKPRqPbs2aOOjg4VFRX1PKc/3EYoSZs3b9Z3\nv/td3Xffffr85z9vRqnAWUPQgu01NDSovr5eEyZM0AcffHDKkBUOh5WUlNTrJZeRSEQ+n08ul6sn\nwMS7+vp6DRw4UNu3bze7lDPSvdZpaWlx/+lut02bNiklJcUy9XYLhUIKhUKWq1uy1o6WdHwnfOvW\nrZaZ7NgtKSlJ77zzjiZMmGB2KaclGo2qqalJe/futdTzQ5KGDBmitWvXWua6he7LjVNTUz92uXH3\nzlY0GtXXv/51JScn9xp+2d2CFVnrlRyIoe5WwZaWFk2dOvW0LiCORqO9Dr1oaWlRbW2tLrjgAvXt\n2/dslHxWNDU1afz48Zb5JFo6/v/dtm3bLLfWjz/+uGbPnm2ptZasO6RBst6OVktLi9ra2ix3ZrK+\nvl4vvviipZ7bx44dU0VFhc4//3xLPUestnsoHV/r8vJy5eTk9Ix2P3FnKxgMauXKlaqurtYf/vAH\nZWZmmlEqEFMELdhWOBxWcnLyKe9VOZ3R7dLx9rtQKKTU1FRVVFScrZLPir///e/q16+fpXYr/H6/\nnE6n5da6rKxMF154oaXWWjp+Pis1NVV79uwxu5QzZrUdre5zWm1tbWaXcka8Xq/Ky8st99wOBALa\ntGmTpQbU9O3bV//4xz80bNgws0s5qY6ODv3Xf/2XmpqaNGTIEK1YsUL9+vWTdPw1MCkpSW63W6+9\n9pqefvppSZLL5eppLezs7NTgwYPldDqVkZGh4cOHS5La2tp04MAB+f1+TZs27SPTgB944AE9+eST\nSkpK0i9+8QtdfvnlkqR169bp9ttv72m5LCkpMXIpAIIW7Kv7hf1UZ6gikYjC4XCvu1jhcLin1z8v\nL88SF6Oe6Le//a1mzZplmU9HDx061LMLF+/n307U3NysuXPnaujQoWaXctqsfJGuZL0dLcl6F0NL\nUk5Ojn74wx9a7nkSjUZVVlamQYMGaciQIWaXc1oqKyv11ltvxe1aL1++XPPnz9f3vvc9Pfzww9q4\ncaNWrlwp6fjP1NraWh0+fFjPPPOM/vnPf8rhcOiiiy7S+vXrlZGRoX//93/XqlWr9Itf/KJnl2ve\nvHmqqKiQ0+nUrbfeqp/85Cc9f195ebnWrFmjsrIyHThwQHPmzFF1dbUk6bbbbtP69es1bNgwFRcX\na/78+ZY504vEYL13hYABuneywuFwr4/x+/3atm2bMjMzVVBQYMmQJVlrtLvP59OePXtUWFhouZB1\n7NgxHT161DJv5rq1t7crPT3d7DJsJT093XIXWmdlZSkUClluYqLD4VBhYaHq6+vl9XrNLue05Ofn\nx/WI95dffrnn8urrrrtOf/3rX3u+5nQ6lZ+fr7KyMk2YMEGpqanKyMjQxRdfrA0bNujgwYM6evSo\nFi5cqL/+9a9KS0vT97//ffn9fo0dO1Zjxoz52N+3du1aLV68WG63Wzk5OcrLy1NpaalKS0uVl5en\n3NxcuVwuLV68WGvXrjVsHQCJHS3gY7qnCp5qdHs4HJbf75fH41FDQ4MaGhoMrjJ2Kioq1NnZGfct\nP9FoVD6fT263W++//77Z5Zyxuro6DR48WO+++67ZpZyRrq4uOZ1Oy73x72a11kHp+HnQpqamk467\njmeDBg3SK6+8ovz8fLNLOWORSESlpaVKS0uL+w9xjh07pqqqKpWWlsZlrY2Njdq3b5/27dun9vZ2\n7d69W2PHjv1IG2F1dbWysrK0efNmud1u1dTU6M9//rPcbreOHj3aM+G0vb1dhw4dUlpamjIzMzVq\n1ChVVVXpy1/+shobGxUIBCRJP//5zyVJDz30kLZu3aolS5YoEolo//79am1tVWZmph5++GFFIhE9\n88wzSk5O7nWoFRBLBC3g/3TvYgUCgV5Ht0vHpxQePHhQkydPttwboRO1trYqEono8ssvj8sf2B+2\na9cuuVwujRgxwuxSPpXDhw+rsLAwbtt9evPee+9p/Pjxlhv9382KrYPBYFA7duzQlClTzC7ljBQV\nFalPnz6We45327dvn/x+vyWCYkpKinJyckwbUDN//nw1NTV97L+vWLFCycnJPc+B5cuXy+VyqaKi\n4iNthO+88478fr+mT5+uTZs2aceOHbrttts0Y8YMfelLX+ppI7zgggs0atQo3XvvvVq0aJGuvPJK\nvf766xoxYoSKiopUUlKi6dOn66WXXtKtt96q//zP/1RlZaXmzp2rrVu36ne/+91HhmpcddVVevLJ\nJw1bJ4CgBej0RrdHo1F1dXUpGo3K4/Fox44dBlcZe+Xl5Tr33HPj/pO9UCikQCCgtLS0k/5wt4I3\n33zTUqPopf8/mMFqo/8/zIo7WtLxgQDxumPRG4/Ho7fffttyd1N9WPel9fF+Pu6cc87Riy++qIkT\nJ5ry9y9fvrzXr/Xr10/r1q1TVlaW/vSnPykjI0NbtmzRiBEj9MADD+jZZ5+Vy+VSfn6+duzYoc2b\nN6t///4KhULyer2KRCKaNGmS3G63Dh8+rL59++qOO+6Q2+3Wvffe2zOgavz48Xruuee0b98+bd68\nWR6PR3379tXEiROVnZ2tBx98UE1NTZo0aZKk4zuB3edjGZIBo8T3KwlggGg0+omj2wOBgHbu3Kmh\nQ4dq2LBhlnrzcyrV1dWaOHFiXH8C3dXVpe3bt2v69OmW3VWRpGeffVZTp06N67U+UVtbmw4dOnTS\ncxFWYcUdLen49+bAgQMtNeL63/7t31RVVWWp5/iJgsGgtm3bpgkTJsT1JMLJkycrJSUlLtf66quv\nVnl5ub73ve+ptbVV3/zmN1VcXNxzPqqiokKrVq3S//7v/yovL0+vv/66jh49qgULFqi5uVnhcFhf\n/OIXdc011+iqq67SI488ooULF+qCCy5Qbm6u9u/fr8bGRt1zzz1atGiRbrzxRv32t79Ve3u7li9f\nrieeeELjx49XVVWVXC6X/vKXvyg7O1t9+/bV888/r1deeaUnnDEkA2cbQQu21d0q6Pf7T/k4n8+n\nsrIyjR49WhkZGYpEIgZVePZ1D8I41dAPs5WXl2v06NFKSkqK6zo/SV1dnS688EJL/Rva2trUv39/\nS9V8IqvWnp6ertbWVksNIhkxYoRee+01y665dHxYQ15ensrLyzVx4sS4/VBt9OjRqqqqiou1vuqq\nqz7SaRAOh9XQ0KDHH39coVCoZ+fo+eefVyQSUXNzs9566y11dnZqzJgx8ng8mjVrlsaOHaunnnpK\n2dnZevXVV7V27VpFo1HdeeedOnjwoFpaWlRfX9/zb7722mv1l7/8Re+++65cLpfGjBmjQ4cOKRAI\naObMmT1Dk6ZPn662tjalp6drwIABWrJkie6++27t379fl1xyifx+v+bMmWOJ7g5YD0ELtpWUlKTk\n5GRt27btlI+LRqNyOBzas2ePJe8ROpWtW7fqoosu+sQ1MFMkElFdXZ3ZZXxmVVVV8vl8cb3WJ4pG\no5Jk6WEv27Zts+SOVvfaW2mKn9frVVVVlaWe472JRCJxPXQnKSlJ7777blys9b333tvr177yla/o\n3XffVVZWllpbWzVkyBDdddddGjdunHbs2KElS5aorKxM0vHv1YqKCnk8Hl1zzTX6wx/+oL59+2re\nvHn6yU9+orS0NDmdTnV0dCgajcrlcunFF1+Uz+eTy+WS3+9XZmamDh48qJqaGrndbiUlJSkjI0Pj\nx4/Xz3/+c+3du1dLly5Vfn5+z92YK1asUHl5uR599FFD1gv2QtCCbSUlJWn27Nlml2GqpUuXasGC\nBZo8ebLZpSS0d2gemgAAFpBJREFUaDSqpqYmLVq0SP379ze7HFu5+OKLLRm0rMjr9epb3/qWZs+e\nzZqfZYMGDdKaNWt08cUXm13KR8yZM0cHDx7s+b3P59N//Md/6Fe/+pUikYiWLFmitWvXauPGjXru\nuee0cuVKzZo1S42NjZo0aZLeeustHTp0SIWFhbrqqqv08ssvq2/fvpo9e7bWr1+vxx9/XKtWrdLk\nyZO1bds23Xvvvbrhhhu0dOnSnj937Nixko5fObB9+/aelsH58+dr69atamlpUUdHh8aPH2/WMsFG\nHN2fmp2mM3owgPgVDofVt29fNTc3q0+fPmaXk9AOHjyoCRMm6PDhw2aXYitJSUkKBAK86TfQueee\nq82bN2v48OFml5LQ/H6/BgwYoGPHjsX14I6WlhYtWrRIe/fu1YEDB7RlyxbNmDFDy5Yt07333quc\nnBzt379fjz32mO6//34dPXpUY8eO1RVXXKGnn35adXV1CgQC8ng8CgaDevPNN7VkyRJ1dHSotbVV\nAwYMUFtbmz744AMVFxfrnHPOUX19vYqKirR3714dOXJEdXV1WrBggbZv3y63263k5GQNHTpUO3fu\nVGZmZs8VFj/60Y909913m71ksI7T6iu25g2rAD6zffv2aeDAgYQsA9TV1SknJ8fsMoCzLicnJyFa\nfeOdx+PR0KFD43at58yZo/Hjx+uiiy7quRMuNTVVDz74oKLRqH72s5/ppptuUnl5uYLBoKZNm6Zd\nu3apoqJC9fX1evvtt5WSkqJBgwbpscce08iRIzVw4EC99tprampq0uTJk/XEE0/I7XZr5MiReuGF\nF+T3+7Vw4ULdfvvtcjgcCofDPWe8rrnmGg0aNEipqalavny5wuGw/vu//1tXXnmlxo0bp7vuukv/\n8z//ozfffNPspUOCIWgBNlVVVWXpaXJWsnv3buXm5ppdBnDW5ebmavfu3WaXYQtjxoxRdXW12WWc\n1IYNG7Rz586P/K+qqkoNDQ3q7OxUJBLRgw8+qJaWFnk8Hp1//vnKy8vTE088oeXLl2v9+vWqqqrS\nihUrtGTJEpWXl8vr9eq+++5TZ2en8vPz9fWvf13hcFhNTU1auXKlIpGI3nnnHd1zzz0qLCxUZ2en\n3G63nnvuOT3zzDNKSUnRuHHjeh67cuVKbdy4UYsWLdLs2bPV3t6uxYsXa9KkSVq5cqXZS4gEQdAC\nbKqqqkoFBQVml2ELdXV1BC3YQm5ubtzusiSaMWPGqKqqyuwyTltWVpZef/11zZs3T3l5ecrMzNRv\nfvMbSdKXvvQllZeXa/Xq1Zo2bZpyc3OVmZmpq666So8//rii0aj+9a9/6ZJLLlFaWpquv/56lZeX\nq6WlRbfddpt++tOfyuVyae7cuUpPT9fLL7+ssWPHas+ePUpOTlZ7e7ueffZZ3XLLLRo1apTq6urU\np08f3XLLLbrnnnsUCATkdDr1s5/9TO+//75WrFhh8mohUcRvYy+As6q6upodLYPs3r1b06dPN7sM\n4KzLycnR66+/bnYZtlBQUKAdO3aYXcYZW7Bgge655x7l5+drwIABGjlypNLT07Vjxw65XC6tXbtW\nSUlJuv3221VcXKwDBw4oPT1dGRkZGj16tN5++23ddNNNamxsVEZGhoqLi+VwOJSVlaXnn39ejz32\nmCT13L3ocDh0+PBhffvb31YgEFBdXZ0cDodcLpd++ctf6oUXXlB7e7uSkpJ06aWXmrw6SDTsaAE2\nReugcWgdhF3QOmiceG4dPJXCwkJNmDBBNTU1uvPOO9WnTx9lZ2frggsu0NKlS7Vp0ybV1dXpkUce\n0c0339zz2pmVlaXhw4fL7XYrPz9f4XBYx44dU1JSkrKzs+V0OjVmzBiFQiGFQiGFw2FlZWXJ5XLp\n0ksvlc/nU0dHh4YPH65gMKiuri7l5+ero6NDycnJCgaDmjFjhubOnaudO3dq6dKlysvL08SJE/Xe\ne++ZvGqwKoIWYFO0DhqH1kHYBa2DxrFa62C34uJi1dTUqK6uTqFQSLW1tZo/f76k4/eXvfnmm/rB\nD36guXPnavXq1WpubtbMmTPlcDg0f/58dXR0qF+/frryyiuVnJysNWvWqLi4WG1tbfrHP/6hhx56\nSG63Wxs2bFAkElFhYaHeeOMNvfrqq5o7d668Xq+WLFmiKVOmqLW1VQ899JDuvvtuLViwQNXV1frO\nd76jyy67TDU1NaqpqdETTzyhb33rWyavGqyK8e6ADXV2dmrgwIE9nwbi7AkEAurXr586Ozvjegxz\nImK8u/EikYjS0tLU2tqqtLQ0s8tJaJFIRP369VNjY6Pl7ud75ZVXdMcdd6izs1NpaWmqqanRihUr\ntGXLFtXW1mrHjh264YYbtG7dOnV2dmrw4MG64447tGzZMrndbg0ZMkTBYFCXXXaZnnvuOQ0fPlx7\n9+6V2+1W//79e5573btYfr9fI0aM0OTJk1VdXa3KykoNGzZMra2tcrvdGjFihKLRqILBoCRpx44d\nmjhxot5//31t3LhRl156qQoLC5WSkqKrr76aM1yQGO8OoDe7du3S6NGjeQNqgPr6emVnZxOyYAtO\np1MjR47Unj17zC4l4TmdTuXn51uyfXDevHmqrq5WfX29wuGw6urq9MMf/lDl5eWaMGGCPB6P1qxZ\no7S0NF1wwQX64he/qNWrV+uRRx5R3759VVpaqo0bN2rjxo1yu9165plnlJSUpKKiItXX18vn8+nw\n4cN64IEH1K9fP7lcLlVVVelzn/ucysvLtXDhQlVUVMjhcMjj8ejYsWOaOXOmNm/erCeeeEIOh0Of\n//zne+rNyMjQU089xaAMnDGCFmBDnM8yDndowW64S8s4Vj2n1S05OVmPPvqoLr/8co0dO1YzZsxQ\nRkaGVqxYoZ/+9KcaPHiwzj//fB05ckT79+/XI488oiFDhkiSxo0bp7S0NHm9Xl177bW65JJLVF5e\nrqSkJM2ZM0d+v1/Lli3T7NmzNWDAAJWWlupb3/qWHA6H/va3v2n06NHq6upSV1eXDh48qEcffVQj\nR47Ut7/9bTkcDs2ePdvk1UEiIGgBNsT5LOMwCAN2w0AM4xQUFFjynNaHde9u1dbWaunSpdq3b59W\nrlyp3NxceTwejRw5Us8++6xuvPFGdXZ2av/+/frxj38sSXK73erXr59qamo0bdo0BQIB5ebm6tVX\nX+25yHjp0qU6evSoFixYoCFDhmjUqFH6/ve/r3379mnKlCk6cuSIurq6tHLlyp6BHNnZ2Vq2bJny\n8/O1bNkytba26qtf/aquuOIKffnLX2ZIBk4bQQuwIUa7G4dBGLAbBmIYx6oDMXrT26CMcDis3/3u\nd/rCF76gBx54QH/6059UXl6u9vb2nkEZgUBAoVBIb7/9tgoKClRfX6+pU6eqT58+OnbsmEpKSlRS\nUqLa2tqetsDJkycrOztbffr00aFDhzR+/HitXr1aQ4cOVSAQUHV1tSZNmqTBgwerrKxM06dP1wsv\nvMCQDJw2Dg0ANlRVVaVvf/vbZpdhC7t379bChQvNLgMwTE5Ojt5++22zy7CFMWPG6OGHHza7jJj5\ncCthZ2ensrKyNG7cON144409Z7BuvfVW/frXv9asWbMUCoX0ta99TZK0efNmpaSkaNasWero6JB0\n/Bzbyy+/rHPOOUc//elPFYlE5Ha71dnZKUlavXq1jhw5IqfTqebmZq1atUpTpkxRVlaWZsyYoby8\nPKWkpMjtdkuSDhw4oLS0NLW0tGjatGk6cuSIGhsbNXToUFPWC/GPHS3AZqLRKK2DBtq9ezdntGAr\ntA4ap6CgQNXV1TrDCdJx7WSDMi677DL5/X7Nnz9fHo9HP/jBD7RkyRLdf//92rBhgySpoqJCc+bM\n0e7du/WPf/xD0WhUBw4cUEVFhYLBoPbs2aP6+nolJSXpvffe07/+9S91dHRo9erVysvL0x/+8Aet\nW7dOX/jCF3To0CH99re/VW1trf7+97+rvb1dklRWVian06msrCxJ0rBhw9TQ0GDaWiH+saMF2ExT\nU5NSUlJ6flDg7Kqrq9Po0aPNLgMwTHfrYDQalcNxWhOQ8Smlp6erb9++amho0LBhw8wuJ6Y+vLvV\n3t6u3NxcjRs3TitWrNDRo0clSTfddJNuuOEG5eXlqaWlRXfddZek44My+vTpo5kzZ+rIkSO65ZZb\neqbsTp8+XatWrVJnZ6dSUlK0fPly1dXV6eabb9auXbtUUlKitWvX9tTx5z//WUePHtV5552nvXv3\n6v777+d5jdPGjhZgM5zPMk5bW5uCwSChFraSnp6ulJQUNTc3m12KLVh98uCpdO9uvfDCCz2voytX\nrtTgwYOVnZ0tj8ejZ599Vrt27dKsWbPkdB5/WxsKheR0OlVXV6c777xTmZmZPX+mw+HQ888/r46O\nDh09elSvvPKKRo8erV//+tfauHGjrrjiCg0ZMkSNjY2SpIULFyovL0/bt2/XokWLNGDAgJ4/a//+\n/crOzjZwRWA1BC3AZmgbNE73IAw+/YTd0D5onEQbiHEyHx6SEQgEtGbNGs2fP/8jj5k/f76eeuop\nSdJzzz2nSy65RA6HQ/Pnz9eaNWvU1dWluro61dTUaOrUqaf8+z78Zz311FO68sore/777373O0Wj\nUf3rX/9Seno657NwSgQtwGa4Q8s43KEFu+IuLeMkwoj3T3LifVuLFi3qaSN88cUXJR1vI2xpaVFe\nXp4efvhhPfjgg5KOtxEuWrRIRUVFuuKKK/TLX/6yp41wyZIlmj59uqqqqjRs2DA9+eSTkqSSkhKt\nX79e+fn52rBhg0pKSiQd32HLzc1VXl6ebrnlFj322GMmrAasxHGGBygT57QlYFPz58/X17/+dS1Y\nsMDsUhLeQw89pMbGxoSaCmYlSUlJCgQCPW+qYJw777xT6enpuvvuu80uJeG99NJLeuyxx/Tqq6+a\nXQpgJ6fVqsKOFmAz7GgZhzu0YFfcpWWcRD6jBVgdQQuwkWAwqPr6eqbgGYTR7rCrnJwczmgZJCcn\nRw0NDerq6jK7FAAnIGgBNlJXV6fs7Oyeyxdxdu3evZsdLdgSwzCMk5KSopEjR2rXrl1mlwLgBAQt\nwEZoGzROOBzW3r17NWrUKLNLAQw3YsQIHThwQMFg0OxSbIH2QSA+EbQAG2G0u3EOHDigzMxMpaam\nml0KYDiXy6VzzjlH+/btM7sUW7DD5EHAighagI1wWbFxaBuE3dE+aBw73KUFWBFBC7ARWgeNwx1a\nsDvu0jIOrYNAfCJoATZC66Bx2NGC3bGjZRxaB4H4RNACbKK9vV3Hjh1Tdna22aXYAndowe64S8s4\nQ4YMUTAYVEtLi9mlAPgQghZgE9XV1SooKJDDcVqXmeMz4g4t2B13aRnH4XBwTguIQwQtwCY4n2Us\nWgdhd7QOGotzWkD8IWgBNsH5LON4vV61tbXp3HPPNbsUwDSDBw+Wz+dTR0eH2aXYAue0gPhD0AJs\ngtHuxtmzZ49Gjhwpp5OXWNiXw+Fg8qCBaB0E4g/vAgCboHXQOAzCAI5jIIZxaB0E4g9BC7CBSCSi\nmpoa5efnm12KLTAIAziOgRjGycvLU21trcLhsNmlAPg/BC3ABhoaGpSenq7+/fubXYotMAgDOI6B\nGMbp06ePBg0apPr6erNLAfB/CFqADdA2aKy6ujp2tACJM1oG45wWEF8IWoANMHHQWOxoAcexo2Ws\ngoICzmkBcYSgBdgAEweNE41GCVrA/8nJydGePXsUiUTMLsUW2NEC4gtBC7ABWgeN09zcLJfLpfT0\ndLNLAUzXp08f9e/fXwcPHjS7FFsgaAHxhaAF2ACtg8ZhNwv4KNoHjUPrIBBfCFpAgvP5fGpsbGQ4\ng0G4Qys+bNq0SdFoVG+//bbZpdged2kZZ8SIEWpublZnZ6fZpQAQQQtIeLW1tcrJyVFycrLZpSS8\n6upqvfDCC/J4PAqFQmaXY2vXXnutotGoFi9ebHYpthYOh+V2u/Xiiy+qsrLS7HISXlJSkvLy8tjV\nAuIEQQtIcLQNGmfdunX605/+pN///vdat26d2eXY2j333CNJKikpMbkSe9uwYYN+85vf6M9//rNe\neukls8uxhYKCAs5pAXGCoAUkOAZhGGfhwoWKRqMqKCjQvHnzzC7H1m688UYVFxfrG9/4html2Npl\nl12moqIiRaNRfelLXzK7HFsYM2YMO1pAnCBoAQmO0e7Gyc7O1pw5c/TMM8/I6eTl1Uwul0ulpaVK\nTU01uxRbczgceuaZZ3TxxRdr1KhRZpdjC0weBOIH7wSABMeOlrHWr1+vSZMmmV0GEDfGjx+vN954\nw+wybIOgBcQPRzQaPZPHn9GDAZgrGo0qKytLlZWVGjx4sNnlAADOstbWVo0aNUrt7e1yOBxmlwMk\nqtP65mJHC0hgLS0tikajGjRokNmlAAAMkJmZKbfbzSXRQBwgaAEJrLttkE81AcA+aB8E4gNBC0hg\njHYHAPthxDsQHwhaQAJj4iAA2A8j3oH4QNACElRlZaVKS0sZqQwANjNq1Cht2bJFFRUVZpcC2BpT\nB4EE9bnPfU4bNmxQcnKyDh48qKysLLNLAgCcZe3t7Ro4cKBCoZAuuugibdy40eySgETE1EHAzhYs\nWCBJuvbaawlZAGAT6enpuuGGGyRJV155pcnVAPbGjhaQoPbu3asLL7xQFRUV6tOnj9nlAAAM4vP5\nVFhYqDfeeEOjR482uxwgEZ3WjhZBCwAAAABOH62DAAAAAGAGghYAAAAAxBhBCwAAAABijKAFAAAA\nADFG0AIAAACAGCNoAQAAAECMEbQAAAAAIMYIWgAAAAAQYwQtAAAAAIgxghYAAAAAxBhBCwAAAABi\njKAFAAAAADFG0AIAAACAGCNoAQAAAECMEbQAAAAAIMYIWgAAAAAQYwQtAAAAAIgxghYAAAAAxBhB\nCwAAAABijKAFAAAAADFG0AIAAACAGCNoAQAAAECMEbQAAAAAIMYIWgAAAAAQYwQtAAAAAIgxghYA\nAAAAxBhBCwAAAABijKAFAAAAADFG0AIAAACAGCNoAQAAAECMEbQAAAAAIMYIWgAAAAAQYwQtAAAA\nAIgxghYAAAAAxBhBCwAAAABijKAFAAAAADFG0AIAAACAGCNoAQAAAECMEbQAAAAAIMYIWgAAAAAQ\nYwQtAAAAAIgxghYAAAAAxBhBCwAAAABijKAFAAAAADFG0AIAAACAGCNoAQAAAECMEbQAAAAAIMYI\nWgAAAAAQYwQtAAAAAIgxghYAAAAAxBhBCwAAAABijKAFAAAAADFG0AIAAACAGCNoAQAAAECMEbQA\nAAAAIMYIWgAAAAAQYwQtAAAAAIgxghYAAAAAxBhBCwAAAABijKAFAAAAADFG0AIAAACAGCNoAQAA\nAECMEbQAAAAAIMYIWgAAAAAQYwQtAAAAAIgxghYAAAAAxBhBCwAAAABijKAFAAAAADFG0AIAAACA\nGCNoAQAAAECMEbQAAAAAIMYIWgAAAAAQYwQtAAAAAIgxghYAAAAAxBhBCwAAAABijKAFAAAAADFG\n0AIAAACAGCNoAQAAAECMEbQAAAAAIMYIWgAAAAAQYwQtAAAAAIgxghYAAAAAxBhBCwAAAABiLPkM\nH+84K1UAAAAAQAJhRwsAAAAAYoygBQAAAAAxRtACAAAAgBgjaAEAAABAjBG0AAAAACDGCFoAAAAA\nEGMELQAAAACIMYIWAAAAAMQYQQsAAAAAYoygBQAAAAAx9v8AXW086dvWKkQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3069ca96a0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "T = stress * u + bound_stress * u_b\n",
    "\n",
    "T2d = np.reshape(T, (g.dim, -1), order='F')\n",
    "u_b2d = np.reshape(u_b, (g.dim, -1), order='F')\n",
    "assert np.allclose(np.abs(u_b2d[bound.is_neu]), np.abs(T2d[bound.is_neu]))\n",
    "\n",
    "T = np.vstack((T2d, np.zeros(g.num_faces)))\n",
    "pp.plot_grid(g, vector_value=T, figsize=(15, 12), alpha=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the traction on face i: T[2*i:2*i+g.dim] is the traction on the face as defined by the normal vectors g.face_normals. This means that for the bottom boundary, the traction T[bot] is the force from to box on the outside (since the normal vectors here are [0,1]), while on the top boundary, the traction T[top] is the force applied to to top faces from the outside (since the normals here point out of the domain)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
